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Abstract 


The project aims at synthesis of a control law for gust alleviation of aircraft. 
The currently favoured Von Karman model of the gust is being used for 
mathematical modelling of the gust. . The gust is assumed to be uniform in the 
spanwise direction and varies only along the vertical direction ie. a vertical 
gust. The statistical properties of gust are described using the power spectral 
approach. The response of the aircraft to the gust thus modelled is studied in 
terms of motions, accelerations, rate of change of accelerations and other 
dynamical characteristics such as passenger comfort. A rigid six degree of 
freedom model of the aircraft is being used whose longitudinal dynamics are being 
considered at present. The performance of the gust alleviation system is evaluated 
by finding the response of the aircraft when it encounters a gust of known 
spectral density specified by the airworthiness requirements. The gust would be 
sensed by sensors such as vanes or alternatively, the aircraft motions sensed by 
accelerometers and gyros and these will be fed as inputs to the gust alleviaton 
system. Either of the methods may be used or their combination. Deflection of 
control surfaces such as elevator, flaps, etc. is being used to counteract the lift 
increments due to the gusts and thus to help in gust alleviation. The deflection of 
controls is governed by the proposed control law ( which is both feedforward and 
feedback) .An exploratory study of the effect of relaxed static stability on gust 
alleviation was also made but the results show that superaugumentation does not 
help in gust alleviation. 

The present project emphasises on the reduction of normal accelerations to 
improve passenger comfort rather than on the reduction of structural stresses. 



The criterion for passenger comfort proposed by Krag is adopted. A combination 
of the feedforward and feedback control law is proposed for the gust alleviation 
system with the objective to improve passenger comfort. This does not deteriorate 
the handling qualities. The feedforward portion senses the gust angle of attack 
and feedforward control is applied through direct lift control devices such as 
flaps while the feedback portion senses the angle of attack and feeds it back to 
the elevator. The best combination of the appropriate gains of feedback and 
feedforward controls is obtained towards minimization of a passenger comfort 
criterion. The determination of the optimal gain values is made by the application 
of DownHill Simplex Method; using the passenger comfort criterion proposed by 
Krag as the objective function. The results show that the combination of 
feedforward loop of gust angle attack to the flaps and a feedback of the angle of 
attack to the elevator reduces the passenger comfort index from 2.17146 for the 
gust alleviation system disengaged to 1.63278 for the gust alleviation system 
engaged which is quite a significant reduction without much affecting the handling 


qualities. 



Chapter 1 . Introduction 


The present work is an attempt for synthesis of a control law for gust alleviation 
of aircraft and to explore the effect of relaxed static stability on gust 
alleviation. 

1.1 Introduction 

The airmass through which airplanes fly is constantly in a state of flux. Its motion 
is variable both in time and space, random and exceedingly complex. This motion is 
called turbulence. The gradients associated with this atmospheric turbulence are 
termed gusts. When an airplane encounters a gust, large sudden loads arise on it 
which may cause failure of the airframe due to accumulation of fatigue, 
deterioration of the handling and riding qualities and especially the 
controllability of the airplane in strong turbulent wind. 

1.2 Literature Survey 

The problem of aircraft response to gusts is one of the the oldest because of its 
strong influence on airplane strength, structural fatigue life, pilot & passenger 
comfort, navigation and aiming, and weapon delivery. The problem was recognized as 
a high priority item by NACA when it was originally established.^ Gust alleviation 
is of considerable interest for small transport aircraft which fly at medium 
altitudes where most turbulence is encountered. For these aircraft, the 
improvement of passenger ride comfort is the main aim of gust alleviation which 
requires reduction of sustained changes in vertical acceleration. 

Because of widespread effects of gust encounter, there has always been a 




continuing interest tor devising means to alleviate or control the loads and the 
motions that result. Essentially, there are two types of gust response control • 
ta; Sensors such as vanes, laser doppler anemometer etc. may be used to sense the 
gusts, its signal being used to control some auxiliary lift de'.'ice i'lhich nullifies 
the forces on the airplane due to the gusts — a typical feedforward process 
(b) some aircraft motion may be sensed by .accelerometers and gyros, which in turn 
IS used to activate some control which counteracts the motion - a typical feedback 
loop process 

The two forms may be used separately or in combination. The direct lift controls < 
EiLC ) such as flaps, spoilers, etc. are more of the gust alleviation type. The force 
applicators (eg. control surfaces) provide either motion control (such as the 
familiar yaw dampers) or may be used to supress modal response. Geometric 
change is limited to very few configurations.^ The structure of the gust response 
problem is conveniently subdivided as is illustrated in Fig i.i.^ 
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Fig i.i. Structure of the gust response problem 


Gust compensation is 


the most efficient method of gust alleviation. A compensation 








of the gust induced lift can be accomplished by a proper deflection of the control 
surface at the very moment the gust begins to act on the wing . At the same time, 
the control surface creates the same lift as the gust but acting in the opposite 
direction,^ Preliminary investigations by Krag^ have shown that for small 
transport aircraft, a mixed open loop and closed loop system might be promising. 
The open loop portion performs a gross alleviation while the closed loop portion 
is responsible for the suppression of residual motions. The gust alleviation 
system affects mainly the short period motion of the aircraft. The vertical 
acceleration is fed back to the control surfaces. Wing control surfaces and the 
elevator are acting in DLC mode. The pitch acceleration generated by the gust is 
not alleviated by this means. This could be accomplished by pitch acceleration 
feedback to the elevator. Both feedbacks are disadvantageous with respect to the 
handling qualities and controllability of the aircraft. This is the reason why © 
feed back is not considered. Because of its static stability, the aircraft tends to 
turn into the direction of the airflow. This behaviour auguments the gust 
alleviation and therefore should not be suppressed by a 9 feedback. An additional 
small 9 feedback would improve the situation without the handling qualities further 
deteriorating. The z feedback has a de-stabilizing effect on the phugoid motion 
and decreases the natural frequency of the short period motion.® According to B. 
Krag, the dynamic response characteristic of the aircraft can be improved by 
feeding back 9 and 9 to the elevator. 

The following methods are also proposed to reduce the vertical accelerations due 
to rough air^ : 

(a) Pitching the whole airplane to maintain a constant angle of attack during 
passage through the gusts. This is advantageous in the sense that it may be 
accomplsihed by the use of the elevator without provision of additional controls 


on the airplane. 



(b) Variation of wing incidence to maintain a constant angle of attack during 
passage through the gusts. 

(c) Operation of flaps or other controls to offset the lift increments on the wing. 

A problem of longitudinal control arises in connection with the methods (fa) & (c). 

A mechanism sensitive to changes in acceleration should be most effective for 
reducing these changes because it could counteract accelerations from any type 
of gust. 

The concept of a power spectrum for the compact representation of a random 
disturbance in time has been successfully used to find gust loads on airplanes^. 
John Tukey has developed a new technique for the derivation of spectra estimates 
from the time history data in which an initial raw estimate of the power spectral 
density is obtained using the Fourier cosine series and this estimate is further 
smoothened by a more effective sharper filter to reduce the aliasing problem 
encountered while finding the power spectrum of the signals in the time domain. 
The author of the same paper also decribes the effect of C.G. position, short 
period damping and wing bending flexibility on gust alleviation. It has been shown 
that the above effects are significant for large flexible airplanes Some of the 
more recent techniques involve an improved mathematical description of an 
aircraft to 3— dimensional isotropic turbulence by establishing a unique 
correspondence between the spatial symmetry properties of the ambient field of 
turbulence and that of the aircraft thereby yielding a more elegant and compact 
formulation offering significant advantages over the direct multiple input 
approach. The new formulation utilizes four describing functions which depend 
upon the gust frequency and the interpanel separation, to determine the normalized 
gust velocity cross spectra for the turbulence field®. ®Eigenspace techniques have 
been used for the design of an active flutter suppression / gust load alleviation 
system for a hypothetical research drone. Full state control laws were designed 



■for the above syatem by selecting feedback gams which place closed loop 
eigenvalues and shape the closed loop eigenvectors so as to stabilize wing flutter 
and reduce gust loads at the wing root while yielding acceptable robustness and 
satisfying the constraints on r m.s. control surface activity. 

Interest in the gust response control has been resumed because airplanes have 
become more sophisticated and versatile in their component buildup that might be 
used to advantage for gust response control purposes without adding complexity or 
compromising the design. Furtherj the need is greater, more is known about the 
problem, and the subject now includes the idea of supressing the structural 
bending effects. Some of the more recent work also incude the concepts of 
controlling the elastic mode response of the airframe called mode stabilization or 
modal response control. It is neccessary for us to continually re-examine and 
improve our analytical and experimental modelling of the atmosphere for 
application to design and flight simulation. Moreover, the recent explosive growth 
in computing capability has provided unprecedented opportunities for 
sophisticated analysis, control and simulation. 

1.3 Active Control Technology 

Active control technology ( ACT ) involves continuous measurement of several 
flight variables via sensors installed at selected stations on the aircraft, complex 
processing of signals obtained from the sensors and feed them back to actuate 
various aerodynamic surfaces to achieve the required objective. ACT is applied in 
some very critical regimes of operation of aircraft involving their safety. Gust 
response and the subsequent alleviation of the loads and motions that result is 
one such area where ACT may be used very effectively because of problems of 
controllability and safety of the aircraft that arise upon gust encounter. 



i.4 Preview of the work 


Chapter 2 describes the characterisation of the gusts. It describes the various 
approaches that have been used so far to model the gusts. The use of the Von 
Karman gust model for the mathematical modelling of the gust for the present 
study has also been justified. 

Chapter 3 discusses how the problem of gust alleviation has been formulated in 
terms of the equations of motion. It also shows the incorporation of a general 
control law — feedforward and feedback into the equations of motion. The signals 
from the feedforward loop are fed to the flap ( direct lift controls ) while the 
signals from the feedback loop are fed back to the elevator. 

Chapter 4 elaborates the method of analysis used in computation and the use of 
power spectral density approach in the analysis. 


Chapter 5 discusses the results obtained from computation. 



Chapter 2. Characterisation Of Gusts 


The influence of turbulence on aircraft flight has always been of concern, since 
the Wright Brothers because of the control upset problem, influence on the static 
design strength of the aircraft and problems in acheiving precise flight. The 
study of gust response of airplanes is a twofold problem requiring the adequate 
representation of the characteristics of atmospheric turbulence and the 
determination of the airplane response ( loads or motions ) in rough air. Trends in 
aeronautics towards higher speed and the development of missiles has served to 
increase the need for more generally applicable techniques, both for the 
measurements of the characteristics of atmospheric gusts and for the calculation 
of the gust response of new airplanes. Our concern is with the dynamic features 
of the problem - the response of the airplane associated with the gradients in 
atmospheric motion which are termed gusts. 

2.1 Mathematical Modelling of Gusts 

For many years, it was common practice to treat the gust response problem as 
response to discrete gusts considered to be isolated encounters with steep 
gradients (horizontal and/or vertical) in the horizontal or the vertical speed of 
air. The earlier gust studies were in terms of the concept of a single or discrete 
gust encounter where the gust shape was considered to be of the sharp edged or 
step function type. The airplane was supposed to be in steady flight in quiescent 
air and suddenly enter a region of updraft. The updraft may be sharp edged or 
graded. Later, this was modified to incorporate a specified profile of the gust. 
Figure Z.lshows the discrete gust approach to design which is still in partial use 


today . 




Fig 2.1. Discrete Gust Approach 


The < i - cosine ) shape of the gust is one of the favoured forms in which Wg is 
the gust velocity and d is the the distance along the flight path. 



Fig 2.2. ( i — cos ) gust shape 

The severity of the gust is controlled by the magnitude of Wm ^nd dm. 

In the light of random turbulence approach, discrete gusts are identified with the 
occasional extreme excursions of the upward velocity. In this sense, the spectral 
approach includes the discrete gusts, while at the same time giving a more 



realistic- arid comprehensive picture. Nowadays, spectral techniques are being used 
for designing aircraft for gust encounter. The assumptions commonly made are : 

(a) Taylor's hypothesis applies, namely that time histories of gust velocities 
obtained from aircraft may be converted to space fixed, temporarily 'frozen', space 
histories 

(b) Gusts are considered to be random in the flight direction only but are 
considered uniform in the spanwise direction due to the notion that the scale of 
turbulence is large as compared to the airplane dimensions. 

The basic spectral procedure involves choosing a spectrum for velocity input, a 
scale value L, an rrns gust severity value cTw and then through the frequency 
response to get the output response spectra as shown in figure 2.3. 
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Fig 2.3. Input Output Relations for gust response 


The frequency response 


function 'H' defined as the response of a variable ' x ' to 



unit sinusoidal gust is fundamental to the spectral approach. Early statistical 
•eatments of the atmospheric turbulence usually employed a spectral shape due to 
'yden for mathematical modelling purposes. This spectrum, based on exponential 


jrve representations of the correlation functions for turbulence velocities. 


chibits an fall off at high frequencies. More recent studies have popularly 

_s 

;ed the spectrum equation due to Von Karman which exhibits an w ® fall off at 
igh frequencies. Figure 2.4 indicates that the Von Karman model of the atmosperic 
jrbulence is the preferred form for modelling purposes and will be used in the 
■'esent study. 



Fig 2.4. Spectral Shape of Atmospheric Turbulence 


'he form adopted by Von Karman for E(fl) is as below 

:(ft) = C (fl/Ap) 2.1 

C 1 + (a/fio> 3 

Phis Exhibits sl bshavioup o*f ft ^ Hnd ft^ low 3.nd hi^h 'fpecfUBnciEs nespsctively. 
rhe constants C and fto and the one-dimensional spectra that are obtained from 


this relation by using the relations in Appendix A are as below ; 






1 

C i + (1.339LA)^]^^® 

ti+|(i.339Lf!.)^] 

C i + (i.339La)^3^^^® 


2.2 


2.3 


Specific items that are usually considered in the description of turbulence are 
the following: 

(a) Proportion P of the total time that the turbulence exists 

(b) Spectral shape ^(fi) 

(c) RMS turbulence severity values z?u> cfy, cj-w 

(d) Integral scale length L 

(e) Peak or level crossing count of turbulence time histories 

(f) Isotropy and homogeneity checks 

(g) Effect of atmospheric variables 

The influence of geographic location, season and altitude on the above is also of 
concern and needs to be considered. 


Turbulence can be further classified into four categories on the basis of 
turbulence severity ; 

( a ) Light 
( b ) Moderate 
( c ) Severe 
( d ) Extreme 

A common means of representing severity of turbulence patches is in terms of the 
rms value ( o- ) of the velocities sensed. The rms value of turbulence, using the 
data on the proportion of time that the patches of various severity are found. 



appears to be essentially invariant with the altitude and is in the range of 3-3.5 
fps. Atmospheric turbulence is a random process in which the maximum values of 
the gust velocities are of the order of 20 times the overall rms value. 

Turbulence severity tends to be linearly related to wind speed U in the form 

tJ’w = c^wo + m U 2.4 

where Cwo is dependent on atmospheric stability ( increasing turbulence with 
decreasing atmospheric stability ) and the slope m appears to be larger for rough 
terrains than for the plains. 

Also intimately associated with the severity and the spectral shape of turbulence 
data is the integral scale of turbulence L, whose mean value lies in the range of 
500—700 ft. Figure 2.5 shows that though the scatter is great, but a rough trend 
indicates that the scale length increases as the severity increases. Press and 
Meadows^® have suggested a value of 1000 ft for L as representative of 
atmospheric turbulence. The worst-case philosophy as specified by the 
airworthiness requirements suggests a value of cr = 10 ft/s. 
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Fig 2.5. Correlation of Integral Scale L with rms gust severity c? 



2.2 Simplifications 


An assumption commonly made is that atmospheric turbulence is isotropic meaning 
that there is no mean rate of transfer of momentum across shearing surfaces or 
Reynolds or eddy shearing stresses, — /3uv, — /Ouw, — jOvw vanish and the cross 
spectra fl^uv; 9 ^uwj ^vw also vanish where u, v, w represent longitudinal, lateral and 
vertical turbulent velocities. Taylor's hypothesis applies, namely that time 
histories of gust velocities obtained from aircraft may be converted to space 
fixed, temporarily 'frozen', space histories. Also, atmospheric turbulence may be 
regarded as a momentarily time frozen random surface in space relative to an 
airplane ie. the properties of stationer ity and homogeneity of atmospheric 
turbulence are also assumed. Statistical functions such as the auto correlation 
function in terms of either time or distance is thus made possible through the 
transformation x = Vt due to Taylor's hypothesis. Some key relations applying to 
isotropic turbulence are listed in Appendix A. 

The overall response for a completely random gust environment is smaller than the 
case when only the random nature of the gust in the flight direction is considered. 
Essentially, the above problem may be treated as that of a system being excited by 
multiple random inputs rather than by a single one. Cross spectra between the 
individual inputs must now be considered along with the spectral value for each 
input. The random velocity of the atmosphere is used as input to the airplane. 
Because of the linearity of the system ( the airplane ), the three component 
velocities of the input can be treated separately. The general gust response 
problem is sub-divided into three subcases : the response of the airplane to (a) 
fore-and-aft gusts Ug, (b) side gusts Vg and (c) up-and-down or vertical gusts Wg. 
Generally, longitudinal gusts may be ignored. Wing loads are primarily governed by 
vertical gusts. Fuselage and tail loads can be obtained by superimposing results 
that are obtained in separate treatments for lateral and vertical gust encounter. 



Each of these may induce technically important airplane response, although it is 
the vertical component Wg that is mainly responsible for normal accelerations 
which IS the main interest of the present study. The analysis of longitudinal 
response to vertical gust requires only a one-dimensional gust spectrum because 
the linear component of the spanwise variation of the vertical gust Wg has no 
effect on the longitudinal motion. Unsymmetrical vertical and side gusts generate 
rolling & yawing moments. They are also generally specified by the airworthiness 


requirements. 



Chapter 3. Formulation Of the Problem 


The treatment of the gust response problem of the aircraft needs a precise 
formulation of the equations of motion. The equations of motion refer to a space 
fixed coordinate system because the gust vel is only defined relative to the 
ground. The longitudinal perturbation equations of motion in steady rectilinear 
flight of a conventional aircraft under the usual assumption of small 
perturbations and linearity ( of perturbed force and moment coefficents as a 
function of perturbed variables of motion ) are well known and may be found in 
standard textbooks on Dynamics of Flight^^'^^ . 


3.1 Equations of motion 

The equations of motion with reference to stability axes are given below in 
dimensional form ; 

u = Xu u + Xj^ u + Xa a + 9 cos Be Q + Xfj fl 3-1 

Ue(a— q) = ZuU + Za a, + + Zqq— g sin 6© ® ^ ^-2 

q = Mu u + u + Ma a + Mj^ a + M^ d + Mq q + }? 3.3 

e = q 

where Ug and ©e airspeed and the pitch angle in the reference flight 

condition respectively; u, a, q, 6 and r? are the perturbation velocity, angle of 
attack, pitch rate, pitch attitude and the elevator angle respectively. The 




dimensional derivatives are defined in Appendix E. 

The following simplifications in the eQuations of motion may be introduced. The 
reference flight condition is a stationary horizontal flight path ($6 = 0). Also 
the thrust T is assumed to be constant so that the derivative Mx and Mx may be 
assumed to be of the order zero. The above equations of motion are non 
dimensionalised and take the form as below. 


( 2MD - 2C|_^tanee - )u - Cx^-a + C|_^e - Cx^T? = 0 3.5 

( 2C|_^- Cz^J )u + ( luD-Cz^B-Cz^ )a - C (2iA +02^)0 -C|_^tanGe 3 6 - = 0 3.6 

"‘“[xiyU — < '“m^s + Cpn^ )a + < igS^ — CmqS )G— = 0 3.7 

q = DO 3.8 


where 




% = 


^yy 

/OSl® 


The non dimensionalized quantities are defined in Appendix B. The relations 
between the non dimensionalized quantities and the longitudinal stability 
derivatives are also given in Appendix B. 

However, for an aircraft provided with pitch attitude feedback to the elevator , the 
equations of motion have to be suitably modified. It may be noted that the feedback 
to the elevator is not the absolute value of the pitch attitude but the 
perturbation value. The above accounts to the addition of effective axial and 
normal force derivatives Cx^ and Cz^ and an effective moment derivative Cm^ to 
the equations of motion. But, since the effect of the elevator on the pitching of 
the aircraft is significantly more predominant in comparison to the axial and the 



normal force derivatives, Cx^ and may be ignored. The term +Cnig S is added to 
the right hand side of the pitching moment equation. 

3.2 Idealization of the airplane 

The airplane may be idealized with varying degrees of approximation corresponding 
to the various assumptions about the way in which it feels the gust structure. The 
simplest approximation is to treat the airplane as a point which flies along a 
straight line. This, evidently, is not a very good approximation. The next step of 
elaboration is to represent the airplane as a segment of x-axis. In this case, the 
variation of the vertical gust in the longitudinal direction is considered but the 
spanwise and the vertical variation is still neglected. This is quite a good 
approximation of the airplane for the present study because normal accelerations 
are mainly due to the vertical component Wg. Reduction of normal accelerations is 
the main criterion to improve passenger comfort. Besides, the analysis of the 
longitudinal response requires only a one-dimensional gust spectrum because the 
linear component of the spanwise variation of Wg has no effect on the longitudinal 
motion. In the case of rigid body motion, no fast motions occur. Also, in a natural ( 
random ) gust field, the low frequency gusts are dominant. Gusts of higher 
frequency appear with small amplitudes as can be seen by the well known Dryden or 
Von Karman gust spectra. In addition, gusts having a smalll wavelength as compared 
with the length of the airplane cannot excite the motion of the airplane because 
they cancel each other. Thus, spectral components of wavelengths larger as 
compared to the length of the airplane l(ie.X^81) are considered where a 
straight line is a fair approximation to most segments of a sine wave. For 
relatively long wavelengths, the variation in Wg along the length of the airplane 
is approximately linear so that the aerodynamic effect of the gust is to modify the 
effective angle of attack and pitch rate. The aerodynamic effect of the gradient 



in iJg equivalent to a certain rate of pitch. By, incorporating aerodynamic force 
terms corresponding to ag and qg into the equations of motion, the effects of 
long wavelength gust components on the airplane response is adequately accounted 
for 


3.3 Modified Equations of Motion 

The equations of motion for controls firmed and a horizontal reference flight path 
after taking the Laplace transform are as follows : 




Wg 

Un 


rrc. )n -4- )r/. — (2Ct4-C~- )0 = s—C- 4-C-^ 




3.10 


(igS‘’~CrnqS)'5 = (~C}*f,^s— Crn^+Cj-nqS 


Wg 

•' uT 


3.11 


where the bar indicates the Laplace transform of the quantities. 

Since direct lift control ( flap ) and elevator control is being used for gust 
alleviation, the control derivatives appear on the right hand side of the above 
equations as : 


= Cx^ ?? + Cxg Sflap 
= T) 4- '^flap 

= Cm^ V + Crr.g Sfiap 


As pointed out earlier in Chapter i, a miKed open loop and closed loop systtsm 



might be promising for small transport aircraft. The open loop portion performs a 
gross alleviation while the closed loop portion is responsible for the suppression 
of residual motions. A general control law was incorporated having an open loop 
control law which presupposes an exact measurement of the gust angle of attack 
signal and a closed loop feedback law which incorporates all the motion variables 
through suitable gains. 


pflap ]ol = ^ i + Ts ) ag 
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where n = g ( q— s a ) 

Thus, r? = ^cl + ^pilot 
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The above is incorporated into the equations of motion as below : 
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where A, B, are defined in Appendix B. 

The above formulation of the equations of motion for the gust response problem is 
used for computation. The ' s ' operator of Laplace transform is written as ik 


where k ( = ftl) is the non dimensional frequency. 



Chapter 4. Method Of Analysis 


The study of gust loads on airplanes is a twofold problem recfuiring the adequate 
representation of the characteristics of atmospheric turbulence and the 
determination of the airplane response ( loads or motions ) in rough air Trends in 
aeronautics towards higher speed and the development of missiles has served to 
increase the need for more generally applicable techniques, both for the 
measurements of the characteristics of atmospheric gusts and for the calculation 
of the gust loads on new airplanes. Turbulence can be defined as the irregular 
random motion of small fluid masses which is so complex that there can be no hope 
of a theory which will describe in detail the velocity and the pressure fields at 
every instant. Existing theories for modelling atmospheric turbulence may be 
classified as either empirical or statistical. 

4.1 Theory of Generalized Harmonic Analysis 

The theory of generalized harmonic analysis® appear adaptable for extending the 
analysis of gust loads beyond the discrete gust case to the case of continuous 
turbulence. Techniques from generalized harmonic analysis involving the power 
spectral density have been used in diverse fields for many years, such as in the ! 
study of random noise problems in communications and in the study of small scale 
turbulence in wind tunnels. Attractive features of spectral analysis for the study 
of gust loads are the possibiities that : 

(a) Continuous turbulence can be described in analytic form by a power spectrum 
rather than by discrete gusts. 

(b) The load response of the airplanes to continuous rough air can be evaluated. 



<c) The desirable response characteristics of an airplane for minimizing gust 
effects in continuous rough air will become amenable to analysis. 

The theory of harmonic analysis is described in Appendix C. 

The response of the airplane to rough air is essentially a statistical one ie. the 
power spectrum of the input ( the air turbulence ) is known and the power spectrum 
of the output ( the airplane response > is required to be found out. From the 
theory of spectral analysis defined above, the link which connects the input and 
the output is the square of the modulus of the transfer functions of the motion 
variables and the gust velocity. The transfer functions obtained have been 
validated by using the numerical example of Sec iO.7 of ^^Etkin. The input output 
relation in terms of the space frequency ( = ^ ) is : 

1 T(iA) 1^ 4.1 

t 

Only the square of the modulus is required, the phase angle of the response is of 
no interest in this connection. The Von Karman model of the gust is being used to 
give us the power spectral density of the input gust spectrum. The longitudinal 
stability derivatives of a Roshkam's business jet transport aircraft given in 
Appendix D are being used for studying the aircraft response to gusts. The data 
for the flap used is also given in Appendix D. 

4.2 Passenger Comfort 

The reduction of accelerations caused by rough air would be of obvious value for 
improving passenger comfort in commercial airline operation. An airplane capable 
of smooth flight through rough air would also be a valuable tool for studying the 
gust structure of the atmosphere. The design of a device for providing 
comfortable flight to the passengers of an airplane in rough air requires a 



knowledge of the factors which contribute to passenger comfort. ^^McFarland 
indicates that slow oscillations of large amplitude are more likely to cause 
sickness than faster oscillations of small amplitude. Some of the results obtained 
in these tests are shown in figure on the next page. 
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The results of these tests showed that very little sickness was produced by the 
shortest period oscillation tested and this appears to be in accord with the 
common experience eg . a periodic motion of a small boat in a rough sea is known to 
produce sickness in a large number of passengers, whereas the motion of a trolley 
car on a rough track, which involves abrupt jolts and jerks containig components 
of oscillation of high frequency and fairly large amplitude, causes sickness in 
very few passengers. Reduction of the more prolonged changes in vertical 
acceleration, thus appears to be beneficial for minimizing sickness. 

Other stimuli that may be important in producing motion sickness include lateral or 
rotational accelerations, motion of objects in the field of vision, and many 
psychological factors such as noise, vibration, temperature, ventilation and so 
forth. Thus, even complete elimination of normal accelerations would not 
neccesarily eliminate airsickness. In flight through rough air, the changes in 
normal acceleration are relatively large, whereas lateral accelerations, rotational 
accelerations, and changes in orientation are relatively small. Peroidic changes in 
normal acceleration are a major cause of motion sickness and its reduction thus 
appears to be the most promising method of improving passenger comfort. Motion 
sickness occurs only at low frequencies ( 0.1 — 0.4 Hz ) in connection with large 
acceleration amplitudes. For the present study, the improvement of passenger ride 
comfort is the main aim of gust alleviation. In order to measure the efficiency of 
the gust alleviation system, the ride comfort criterion is expressed in terms of 
the motion parameters. This can be accomplished by definining a comfort criterion. 
The main parameter influencing passenger comfort is the vertical acceleration. 
The root mean square ( r.m.s. ) value of the vertical acceleration is given by 

C = z = 4^ 


4.2 



In the time domain, 


T 

= lim ^ (t) dt 

T -too ^ I J 

-T 

where T is the time of observation. 

In the frequency domain, 



(cj) dw 
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where *^he power spectral density of the vertical acceleration. 

The above equations represents the simplest of the comfort criteria. In terms of 
the above index of passenger comfort, reduction of the root mean square ( r.m.s. ) 
values of the vertical acceleration at various locations of the aircraft cabin 
would improve passenger comfort. 


4.3 Computational Procedure 

The study of gust response of aircraft in longitudinal motion has been done 
computationally as follows : 

( i ) Computation of the roots of the characteristic equation. 

( ii ) Computation of the time responses of the motion variables like u, a, q, 6 and 
load factor n. 

( iii ) Computation of the frequency response ( power spectra ) of the motion 
variables through the power spectral density approach. 

( iv ) Reduction of the r.m.s. value of normal accelerations through the proposed 
control law by optirnising the control law variables ( DownHill Simplex Method ) to 

improve passenger ride comfort. 

( V ) Includes the idea of a random gust by simulation of a random gust through a 



control law by optimising the control law variables < DownHill Simplex Method^® ) to 
improve passenger ride comfort. 

( V ) Study the effect of relaxation of static stability on the gust alleviation of 
the aircraft . 



Chapter 5. Discussion Of Results 


This chapter discusses the results of the computations presented in the tables 

and figures. The effect of various feedforward and feedback laws is discussed. An 

4 

exploratory study into the effect of relaxation of static stability and its 
artificial compensation or superaug. mentation through pitch attitude feedback to 
the elevator is discussed. A comparative study of the various control laws 
obtained by the minimization of the objective function ( passenger comfort index ) 
w.r.t. the various gains for the controls to be applied, is made through Table 5.1. 


Table 5.1. Comparative Study Of The Various Control Laws 


Control Law 
Gains 

1 Gust alleviation 
system disengaged 

2 Kccg = 5.250667 
'^filter = 358.006165 

3 Kq = 0.721041 

4 Ku = 4.503037 
Ka = 0.697434 


5 K 


ag 


17.842354 


Tf liter = 83.386169 
Ka = 0.134872 
6 Kcc = 0.057423 


Comfort index 


2.17146 


2.03345 

1.77057 

1.65529 

1.63278 


Short Period 

Frequency ( Hz ) thaif ( sec ) 
2.65 0.685 


2.65 

4.3573 

4.4185 

3.065 


0.685 

0.842 

0.633 

0.682 


2.05263 


2.8341 


0.684 



Fig 5.i shows the power spectrum of the atmospheric turbulence as obtained from 

the Von Karman model of the gust as a function of the space frequency a = ^ 

Ve 

Fig 5.2 compares tht; various gust alleviation systems in terms of the reduction of 
load factor. Since the gust alleviation system affects only the short period motion 
of the airplane , we shall consider reduction of the load factor near the short 
period frequency. Refer Table 5.1 for the explanation below. 

CASE 1 ; When the gust alleviation system is disengaged, the r.m.s. value of the 
load factor ( a measure of passenger comfort ) is 2.17146. The conventional 
airplane has been assumed to have good handling qualities. 

CASE 2 ; When the gust alleviation system 2 ( a feedforward loop which feeds the 
gust angle of attack ag to the flap after passing it through a proportional plus 
derivative filter ) is engaged, the r.m.s. value of the load factor reduces to 
2.03345 ( ie. 6.5'/, ) which is not very significant. The frequency and the time to 
halve the oscillations of the short period remain the same. This means that there 
is no change in the handling qualities because damping of the short period mode is 
a crude estimate of the handling qualities. 

CASE 3 : When the gust alleviation system 3 ( pitch attitude 0 feedback to the 
elevator ) is engaged, there is a greater reduction in the r.m.s. value of the load 
factor than Case 2, a reduction of 18.6% which is quite significant with reference 
to our main aim of improving passenger comfort or alternatively, decreasing the 
r.m.s. value of the load factor. The frequency and of the short period mode 

increases to 4.3573 Hz and 0.842 sec respectively. The increase in the short 
period frequency is helpful for improving passenger comfort because low 



frequency oscillations are mostly responsible for causing motion sickness in most 
passengers. The decrease in short period damping indicates a deterioration in 
handling qualities. 

CASE 4 : When the gust alleviation system 4 ( forward velocity u and angle of 
attack a feedback to the elevator ) is engaged, a reduction of 23.8S'/. in the r.m.s. 
value of the load factor is achieved. The frequency of the short period mode 
increases to 4.4185 Hz while decreases to 0.633 sec. The increase of the 

short period frequency is helpful for improving passenger comfort. The increase 
in damping of the short period mode is not very much and therefore not much 
deterioration in the handling qualities. Thus, this gust alleviation system appears 
to be promising towards improvement of passenger comfort. 

CASE 5 -■ When the gust alleviation system 5 ( a combination of the feedforward 
and feedback loop — the gust angle of attack ag is directly fed to the flap after 
passing the signal through a proportional plus derivative filter and feedback of a 
to the elevator ) is engaged, a reduction of 24.91% in the r.m.s. value of the load 
factor is achieved. The frequency of the short period mode increases to 3.065 Hz 
while t^alf decreases to 0.682 sec. The shift of the short period towards the high 
frequency regime helps in improving passenger comfort. The increase in short 
period damping is very negligible and therefore no deterioration in the handling 
qualities. From the above, this gust alleviation system is more promising towards 
improvement of passenger comfort without the deterioration of handling qualities. 

CASE 6 : When the gust alleviation system 6 ( angle of attack a feedback to the 
elevator ) is engaged, the r.m.s. value of the load factor reduces to 2.05^63, a 
reduction of 5.6% which is an insignificant improvement in passenger comfort as 



compared to Case 5. The frequency and of the short period mode increases to 

2.8341 Hz and 0.684 sec respectively. 

Fig 5.3 shows that significant changes in forward speed are all associated with 
the phugoid mode as expected since speed changes are very small in the short 
period mode. 

Fig 5.4 shows that the a response to gusts is flat except for a rise at the short 
period after which it falls of rapidly. The a response of the gust alleviation 
system 5 shows that the response is lesser than when the gust alleviation system 
is disengaged. 

Fig 5.5 shows the pitch rate response. The response of the pitch rate becomes 
lesser when the gust alleviation system 5 is engaged than the case when gust 
alleviation system 4 is engaged or the system disengaged. 

The load factor response is plotted in Fig 5.6 The load factor reduces when the 
gust alleviation system 5 is engaged and the r.m.s. value of the load factor 
reduces by almost 25 ’/, which significantly inproves the comfort of the passengers. 

Fig 5.7 to Fig 5.10 show the output spectra obtained by multiplying the 
corresponding plots by the Von Karman gust input of Fig 5.1. Their shapes are 
similar to those of the reponse functions except that they fall off more rapidly at 
the high frequency end. 


Fig 5.11 shows the time history of a random vertical gust simulated through random 
number generation^*^”^^ to have the power spectrum similar to the Von Karrrian gust. 



"ig 5.12 ihows the power spectrum of the random gust and compares it with the Von 
Carman gust model. 

"ig 5.13 to 5.16 show the output time history of the motion variables for the 
■'andom gust for the cases when the gust alleviation sytem is engaged and 
disengaged. 


-ig 5.17 shows the gust alleviation system. 
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Fig 5.12 Power Spectrum of the Simulated Random Gust 
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Fig 5.1-^ Time History of <(1) for the random gust 



Fig 5.15 Time History of ^(t) for the random gust 



Fig 5.16 Time History of Tt (t) for the random gust 
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Fig S']7 Gust Alleviation System 







Conclusions 


The improvement of passenger comfort is the main aim for gust alleviation of 
aircraft for small transport airplanes. The efficiency of the gust alleviation 
system is measured in terms of a simple comfort index adopted by Krag^. In view of 
the above, it is concluded that 

( i ) A combination of an open loop ( incorporating direct measurement of the gust 
angle of attack and a proportional plus derivative filter to the direct lift control 
(flaps) device) and a closed loop ( incorporating a feedback of the angle of attack 
a to the elevator ) is promising for the gust alleviation of aircraft. The 
reduction in the comfort index is quite significant, more than other cases without 
much change in the handling qualities of the aircraft. 

( ii ) A closed loop ( incorporating a feedback of u and angle of attack a to the 
elevator ) also appears promising w.r.t. gust alleviation of aircraft as the 
reduction in the comfort index is significant but there is a deterioration in the 
handling qualities of the aircraft due to significant increase in the short period 
frequency. 

( iii ) An exploratory study into the effect of relaxation of static stability and its 
artificial compensation by pitch attitude feedback to the elevator does not help 
in alleviation of loads and motions that result, when an airplane encounters a 
gust. 

Thus, a combination of feedforward and feedback loop process is helpful for 
gust alleviation without much change ( deterioration ) in the handling qualities of 


the airplane. 



Suggestions For Further Improvements 


'hs present work on the ntudy of guzst alleviation of aircraft rriay be extended as 
jelow : 

1. The study of gust alleviation of aircraft may be extended to both 
symmetric and asymmetric rnaneuvres. 

2. The effect of time delay which is inherent in complex digital control 
systems may be accounted for and studied. 

3 . The study would be more realistic by the inclusion, of 3-d turbulence and 
,aking into account, the complete set of the equations of motion. 

.4. A more comprehensive index of passenger comfort should be used to 
iCGQunt for various factors such as noise, location etc. 


f*' XV J 
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Appendix A 


( i ) Isotropic Turbulence relations 

Some key relations relating to isotropic turbulence •• 

1-D spectra 
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The i-D spectrum <^.u is defined in terms of a general 3-D spectrum E(Jl); fi= 
where to = circular frequency and V = speed of flight. The vertical spectrum is 
expressible in terms of the longitudinal spectrum ^u- 

The auto correlation function R^j and the spectrum are the Fourier transform 
pairs but because they are both symmetrical can be expressed simply by cosine 
transforms. The area under the spectrum represents the mean square value of 
turbulence. The initial value of the auto correlation function by definition is also 
the mean square value. The integral scale L is defined as twice the area under the 
correlation function normalized by cr^^. The initial value of the spectrum ( 
using the frequency argument ^ ) is 
The form adopted by Von Karman for E(iT) is as below 


E(A) = C 


C i + 


A.ii 


zi 

This exhibits a behaviour of A"* and A® at low and high frequencies respectively. 
The constants C and Ap and the one-dimensional spectra that are obtained from 
this relation by using the above relations are as below : 
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( ii ) Simulation of random gust 

zo-2iThe random gust is simulated by considering the vertical gust as a Gaussian 
random process Wgtt) with mean zero and the mean square spectral density function 
E'o^w). This process can be simulated by way of the following series : 
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(? is defined as the standard deviation of the process Wg(t), Uf, < k = 1, 2, N), 

are independent random variables identically distributed with the density function 
g(a)) = g(U[^) obtained by normalizing So(c«)) ; 


cj(w) = §2^ ; A.ie 

CT*" 

and are independendent random variables identically distributed with the 
uniform density ^ 5 ^ between 0 and 2 it. Note that and ( k/1 = i, 2, N ) are 


independent . 
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Appendix B 

( Non Dimensionalization ] 


The dimensional derivatives used in the equations of motion are defined as below: 

Table E.i. 
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Non dimensionalized stability derivative relations 
The various non dimensional systems in use vary in : 

( i ) The definition of the time constant t 

( ii ) The choice of the characteristic length 1 

( iii ) The definitions of the non dimensional stability derivatives 

The system followed throughout the oraseht Hork uses the mCA sUbility 

derivelives »ith time constant I* = i and the characteristic length used in the 

longitudinal equations of motion as 1 = f ■ The relative mass or density parameter 

U defined as the ratio of the airplane mass to the volume SI of air is expressed 

^ ~ Mi ' 

The non dimensional system used is defined in Table B.2 



Table B.2. 
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t 
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Jl- 

D =. -i 

dt 


dt 


The non dimensional derivatives are related to the stability derivatives given for 
the aircraft as in Table B.3 below : 


Table B.3. 
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Czrj = 
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The nomenclature used in the equations of motion after incorporating the suitable 
modifications and used in the computer program is as below 



Table B.4. 
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E3 ; 





B = -Cxu - ZCL^tanOe - Cx^Ku 

^ — Cxj^Kcc 

F = Cl^ - Cx^Kq 
H = 2/A - Cz^ + Cz^KnT^ 

^ ~ (2/A+Czq) — Cz^^ (Kq + KnT|^) 

L = — Cpfiy — Cm^Kq 
0 = — — Cm^Ka 

Q ~ l^nT^) 

V = Czq - Cz^ + Cz^KagT-fiiter 
X = Cmq - Cm^ + 


T = CxjKagTfiiter 



Appendix C 

( Harmonic Analysis Theory ) 


The theory of harmonic analysis indicates that an arbitrary periodic function can 
be represented by a Fourier series. In order to apply this technique to non- 
periodic phenomena, the Fourier series takes the form of the Fourier integral. A 
neccessary condition for the frequency representation Gtu) of a time function F(t) 
to exist, is that the Fourier integrals must be convergent. In order to develop a 
frequency representation applicable to continuous disturbances, the theory of 
random processes makes use of the concept of a stationary random process ie. the 
disturbance is invariant with time and a statistical equilibrium exists. For a 
stationary random function y(t), the mean square y^(t) represents a measure of the 

disturbance intensity. 




C.i 


The funolion a<l) is considered to be ocmoosed of an infinite number of sinuso^ 
components with circular freauenca o, between 0 and The portion of a (t) 
arising from the oomoonents hading frequenoies between a and a+du is denoted by 
^(u)du. The function #(u) is called the power spectral density f p. s. d. ) function. 
Atmospheric turbulence is a stationary random process that can be statistically 

described by a single reduced p.s.d. curve. 


oo 


y®(t) = J ^(w)du 


{ 


dtr 


C.2 


C.3 



where I 1 indicates the modulus of the complex quantity. Thus 4>iLi)dw is a measure 
of the contribution of that frequency to y'^(t). The system frequency response T(iu) 
is defined such that T(iu)B^‘^^ is the system response to the sinusoidal input 

oo 

T(iw) = iu J A(t) e"^'^ dt 
0 

where A(t) is the response to unit step. The p.s.d. of may be found by first 

determining the airplane response to a step gust and then using this step 
response in equation (8) above to determine the frequency response for a 
continuous sinusoidal gust input. 

The auto correlation function R(t) is defined as 


R(t) = Lim i f y(t) y(t-t-T) dt C.5 

T -400 1 J 
0 

and it is reciprocally related to the p.s.d. function as 

OO 

r(t) = J cos(wt) du C.6 

0 

00 

== I J R(t) cos(wt) dr 
0 

A change of variables is appropriate in the gust case because in these terms 
power spectra are independent of the airplane forward speed. Use is made of the 
distance x in feet and the associated frequency arguments Si. Thus the turbulent 
vertical velocity distribution along a line in space can, at an instant of time, be 
considered to represent a stationary random function of space consisting of an 
infinite number of harmonics of various frequencies or wavelengths. 



X = Vt 


C.8 


Ths vsiri3,bls is ths rsducsd 'fpsciusncy in psdisns psp “Foot snd vHPisbls x is ths 
aipplane flight distance in feet. Thus, 

= V^(VA) Q g 


The resulting power spectrum of atmospheric turbulence is independent of the 
aipplane speed and in a limited sense may represent a universal reduced power 
spectral density function of atmospheric turbulence. 

oo . 

t(ifi) = in J A(x) dx C.iO 

0 


A 

where A(x) = A(t> = response to unit step gust input 
The input output relation in terms of n ca be written as 

I T(iA) 1^ C.ii 

The probability distribution of the output can be represented by a normal 
distribution with mean 0 and the standard deviation o* obtained as below 

QO QO r 

= J ^o<w> dw = J ^o<n) dn 

0 0 

This relation ties the normally distributed loads to the time history and is thus of 
importance in the gust load analysis 



Appendix D 

( Longitudinal Derivatives of Example Aircraft ) 


Longitudinal Aerodynamic 
Aircraft 
W = 13000 lbs 

lyy = 18800 slug ft^ 

Xcg = 0.315 


Characteristics 

h = 40000 ft 
M = 0.7 

S = 232 ft^ 


of Business Jet Transport 

i 0 = 0.000588 slug ft'® 
c = 7.04 ft 

00 = 0 (stability axes) 


CC 0 = 2 .7 


0 


Non Dimensional Cruise Non Dimensional Cruise 


Derivative 

Ue= 675ftsec'^ 

Derivative 

Ue = 675 ft sec' 

^Du 

0.0 

C| 

0.556 

^De 

0.033 

Cmu 

0.05 


0.0 

Cme 

0.007 

Ct 

' Xe 

0.033 

^"’Tu - 

0.0034 


0.30 


0.007 

^Le 

0.410 

^mct 

0.64 


0.0 

Cmj 

0.0 

^Lu 

0.40 


6.7- 


5.84 

Cmq 

15^0 


2.20 


L32 

C| 

4,70 



For the flap used : 

Cd = 0 

^Sf 

^Ljfr ~ ^Lq 

“-Sf ® 





Appendix E 
( Program Listing ) 


'ogram thesisdnput, output), 
:Qnst 


elevatDr_deflection=0.0; 
delta_tirne=D .074; 
lnitial_time=0 .0; 
Final_time=9.0; 
lnitial_ Altitude=40G00 .0; 
Initial_omega= i E-5; 

F inal_omega= iE-2; 

gr=32.2; 

pi=3. 14 15926; 

np=7; 

mp=8; 

ftol=iE-6; 

nnk=3; 


nn=256; 

nnby2=128; 

nn2=512; 

0rder_m=4; 

Isign=-1; 


eps=lE-6; 

Total_N=256; 

type 

perturbation_quantity = (u,w,q, theta); 

matrix = array [perturbation_quantity] of real; 

Planeinf oT ype=(Ve,Rho, 

S,c,Iyy,thetaE,Muul, 

CLe,CLu,CLalpha,CLalphaDot,CLq,CLdeltaE, 

CDe,Cte,CDu,CTu,CDalpha,CDalphaDot,CDq,CDdeltaE, 

CMu,CMalpha,CMalphaDot,CMq,CMtheta,CMdeltaE); 

DerArr = arrayiPlanelnfoTypel of real; 

glmpnp = ARRAY [l..mp,l..np] OF real; 

gimp = ARRAY Cl..mp] OF real; 

glnp = ARRAY [l..np3 OF real; 

double = real; 

galr = array[1..97] of real; gldarray = arrayCi..nn2] of real, 



glnarray = arr3.yCi..Tolal_N] of real; 

glmarray = array[i..Order_m] of real; 

var 


PlaneProg 

■Der Arr; 

CXu,CX,CXalpha,CZu,CZjCZ3lpha,CXq,Czq real; 

CZalphaDot,CXalphaDot,CXd£ 

iltaE,CZdaltaEreal; 

CXdeltaFlap,CZdeltaFlap 

real; 

CMdeltaFlap 

real; 

DskFLpil,fli;kii,inp 

■text; 

kki,kk2,kk3,kk4,kk5 

:text; 

p 3 d_w,big_omega 

real; 

Frequency, dim_frequency 

real; 

Sigma_w,Scale_turbulencB 

real; 

response 

■matrix; 

time 

real; 

FWg,Wg,Wmax 

real; 

Load_factor 

■real; 

Ph_u_gust,Ph_w_gust,Ph_q_ 

.gust real; 

Ph_theta_gust,Ph_LF_gust 

real; 

Ph_u_elev,Ph_w_elev,Fh_q_ 

elev real; 

Ph_theta_elev,Ph_LF_elev 

real; 

TF_u_gu5t,TF_w_gust,TF_q. 

.gust real; 

TF_theta_gust,TF_LF_gust 

real; 

TF_u_elev,TF_w_elev,TF_q_ 

.elev real; 

TF_theta_elev,TF_LF_elev 

real; 

Dut_w_gust,Out_q_gust 

real; 

Out_LF_gust 

■real; 

Out_u_gust,Out_theta_gust 

real; 

ACMu,ACMalpha,ACMaACMdeltaE real; 

ACMalphadot,MuUjt_star,iB 

real; 

ACMthela 

■real; 

Ch_length,ACLe,CLetanthetaE real; 

Ku,Kalpha,Kq,KthetaXLdFctr 

real; 

Kalpha_gust,T_F liter 

real; 

Aa,Ab,Ac,Ad,Ae,Af,Ag,Ah 

real; 

Ai,Aj,Ak,Al,Am,Ao,Ap,Aq 

real; 



Ar,Gt,Gu,Gv,GWjGx,Gy 

real; 


E1,E2,E3 

real; 


TFactor.KF actor_CLaw 

real; 


sumi,surn2,corrifort_criterion .real; 


delta_orriega 

real; 


check_coeff, count 

real; 


temporarily 

real; 


i,iter,indim 

Tnteger; 


X 

•glnp; 


y 

■■gimp; 


p 

:glmpnp; 


pr 

:3lno; 


□rnegak;Phik,T ime_Data 

:glnarray; 


random_Phi,random_omega 

■■real; 


random_T ime_F unc 

real; 


accept 

:boolean; 


Iomega, Iphi,lrandom_y 

integer; 


iicount,nkki 

:integer; 


gomega,gphi,grandom_y 

:galr; 


golixi,golix2,golix3 

integer; 


gplix i,gplix2,gRlix3 

integer; 


gylixi,gylix2,gylix3 

integer; Moutu,Moutw,Moutq,MoutLF 

:gldarray; 

wki,wk2 

:glnarray; Mouttheta 

■gldarray; 

wkm 

;glmarray; Frequencys 

real; 

fdt,pm 

real; temp 

integer; 

cof 

:glmarray; 


nnkl,kkm 

integer; 



FUNCTION sngl(xx:reaD:real; 
BEGIN 
sngl := XX 
ENDj 

Procedure ReadDekForAcDer; 
var 


1 


:PlaneInfoTypej 



Procedure ReadDskForAcDer; 


var 

i ;PlaneInfoType; 

begin { ReadDskForAcDer } 
for i:=Ve to CMdeltaE do 
begin 

Readln(DskFl,PlaneProgCi]); 

end; 

end; C ReadDskForAcDer } 

procedure Am_ConstCoeff; 
begin C Am_ConstCoeff} 

CXu ;=PlaneProgCCT u]-PlaneProg[CDu]; 
CX:=PlaneProg[CTe3-PlaneProgCCDe]; 
CXalpha:=PlaneProg[CLe3-PlaneProg[CDalpha]; 
CZu ;=-PlanePr ogCCLu]; 

CZ:=-PlanBPro 9 CCLe]; 

CZalpha:=:-(PlaneProg[CDe3+PlanBProg[CLalpha]); 
ACMu :=PlaneProgCCMu]; 

ACMalpha ;sPianBProgECMalpha3; 

CXa ;=-PlaneProgCCDq3; 

CZq:=-PlaneProgCCLa3; 

ACMq -'sPlaneProgECMq]; 
CZalphaDotte-PlaneProgECLalphaDot]; 

CXdeltaE :=-PlanePragCCDdeltaE3; 

CZdeltaE ;st-PlaneProg[CLdeltaE]; 

ACMdellaE ;=PlaneProgCCMdellaE3; 
CXdeltaFlap:=0.0; 

CZdellaFlap ;=-0 .2*PlaneProg[Cle]; 
CMdeltaFlap;=0.0; 

CXalphaDot:=-PlaneProgCCDalphaDot]; 

ACMalphadol:=PlaneProgtCMalphaDot]; 

ACMtheta :=PlaneProg[CMtheta3; 

Muu :saPlaneProg[Muul3; 

l_slar:=PlaneProg[c3/(2fPlaneProgIVe3); 

Ch_length :=PlaneProgCc3/2; 



iB:=PlaneProgClyy]/(PlaneProgCRho]fPianeProgCS3* 

sqr(Ch_length)*Ch_length); 

ACLe ;=PlaneProgCCLe]; 

CLetanthetaE •=PlaneProg[CLe]*sin(PlaneProg[lhetaE]) / 

cos(PlaneProgCthetaE3); 
end; { Am_ConstCoeff 3 

procedure Dummy_Coeffs; 
begin { Dummy_Coeffs } 

TF actor :=C2*sqr(PlaneProgC Ve]))/(grf PlaneProgCc]); 

Aa:=2*Muu; 

Ab:=-CXu-2*CLetanthetaE-CXdeltaEfKu; 

Ac :=:CXdeltaE*TFactor*KLdFctr; 

Ad :=-CXalpha-CXdeltaE*Kalpha; 
Ae;=-CXdeltaE*(Kq+TFactor*KLdFctr); 

Af :=ACLe-CXdeltaE*Kthela; 

Ag :=2*ACle-CZu-C2deltaE*Ku; 

Ah '-= 2 ^ Muu-C2alphadot+CZdeltaE^KTFactorf KLdF ctr; 

Ai :=-CZalpha-CZdeltaE5KKalpha; 

Aj:=-(2fMuu+CZq)-CZdeltaE*(Kq+TFactorIKLdFctr); 

Ak:=-CLetanthetaE-CZdeltaE*Ktheta; 

Al :=-ACMu-ACMdeltaEf Ku; 

Am:=-ACMalphadot+ACMdeltaE*TFactor5KKLdFctr; 

Ac ACMalpha- ACMdeltaEf Kalpha; 

Ap:=iB; 

Aq:=-ACMq-ACMdeltaE*(Kq+TFactor*KLdFctr); 

Ar :=:-ACMtheta-ACMdeltaEf Ktheta; 

Gt:=CXdeltaFlapfT_Filter*Kalpha_gust; 

Gu:=:-CXalpha+CXdeltaFlap*KalDha_gust; 

Gv:=CZq-CZalphadot+CZdeltaFlapfT_Filter*Kalpha_gu5t; 

Gw:=-CZalpha+CZdeltaFlap*Kalpha_gust; 

Gx := ACMq- ACMalphadot+CMdeltaFlap*T_F ilter*Kalpha_gust; 
Gy:=-ACMalpha+CMdeltaFlap*Kalpha_gust; 

Ei:=CXdeltaE; 

E2:=CZdeltaE; 

E3 *=ACMdeltaE; 



end; C Dumrny_Coeffs } 


procedure F errari(A,b,c,d,e real; 
var rei,re2re3re4,iml,im2,im3,im4 real); 
var 

Pjq,R,S,smr,sms,mu real; 

Procedure Quad(a,breal; var Re,Im,Rel,lmireal); 
var 
Dreal; 

begin { Quad } 

D;=a*a-4fb; 
if D<0 then 
begin 

D:=sqrt(-D); 

Re;=-a/2; 

Rei:=Re; 

Im;=D/2 ; 

end 

else 

begin 

D:=sqrt(D); 

Re~( -a+D)/2; 

Rei.'=( -a-D)/2; 
lm;=0; 

end; 

end; C Quad 3 
Procedure NR(var xreal); 

Function f(x real) real; 
begin 

f ••=:xlx*x-2f cf sar(x)+Sf x-R; 

end; | 

Function fprime( var xreal) real; 
begin 

f prime :=3fx*x-4Ic*x-+'S; 

I 



end; 

begin 

repeat 

X :=x-F (x)/F prime(x); 
writeln('inside NR'); 
until abs(F (x)/F prime(x))<0 .0000001; 
end; 

Procedure Bab( var mu real ); 
var 

VjOldMureal; 

begin 

repeat 

writelnC'inside Bab'); 

01dMu;=mu; 

V :=S+01dMu*(01dMu-2*C); 

Mu:=R/v; 

until abs( mu-OldMu )< 0.0000001 ; 
end; 

procedure Birge_Vieta(var xreal); 
var 

aal :arrayCi..nnk3 of real; 
bbi,ccl :arrayC0..nnl<] of real; 
kk integer; 
begin C Birge_Vieta 3 
aal[13:=-2*c; 
aalC2]:=S; 
aai[33:=:-R; 

while (bblCnnk3 > eps) do 
begin 

writelnC'inside Birge_Vieta'); 

bblC03:=l; 

cci[03.=l; 

for kk:=l to nnk do 
begin 

bb iCkk3 :=aa lCkk3+x*bb ICkk- 13; 
cc lEkkl :=ibb l[kk3+x*cclCkk- 13; 



end; 

x--x-bbiCnnk]/ccl[nnk-l]; 

end; 

end; C Birge_Vieta ) 
begin C Ferrari } 

S:=b*d+c*c-4*afe; 

R;=BfC*D-B*B^^E-AfD*D; 

mu:=R/S; 

CNR(mu); 

Bab(mu);} 

Birge_Vieta(mu); 
if not( mu< sqr(b)/(4^l^a) ) then 
writelnCf 11/The roots of oh eqn do not exist') 
else 
begin 

P:=( (B/A) +sqrt( sqr(b/a)-4)Kmu/a) )I0.5; 
q;=(b/a)-p; 

smr;=( p*((c/a)-p*q)-(d/a) )/(p-q); 
sms:=c/a - p)(Eq-smr; 

Quad(p,smr,re i,imi,re2,im2); 
Gluad(q,sms,re3,ifn3,re4,im4); 
end; 

end; £ Ferrari } 

procedure New_ChEqn(var Al,Bl,Ci,Di,Ei:real); 
begin £ New_ChEqn } 

Ai:=Aa«Ap*Ah; 

Bi.-=AbfAp*Ah+Aa*(Ah*Aq+ApittAi-Am«Aj)-Ag^)fAc*Ap; 
Ci:=Ab*(Ah*Aq+Ap*Ai-Am*Aj)+Aa*(Ah*Ar+Ai*Aq-Am* 
Ak-Ao)lcA3)-Agf<Acf Aq+AdfAp-Am*Ae)+Al*(Ac*Aj- 

Ah*Ae); 

Di:=Ab*(Ah*Ar+AilAq-Am*Ak-Ao*Aj)+Aa*(Ai5tAr-Ao* 

Ak)-Ag*(Ac*Ar+Ad*Aq-Am*Af-Ao*Ae)+Al*(Ac*Ak+ 

Ad)fcAj-Ah*Af-Ae*Ai); 

Ei:=Ab*(Ai*Ar-Ao*Ak)-Ag*(Ad*Ar-Ao*Af>+AH(AdIAk- 

Ai*Af); 



writeln(fii/The characteristic eqn is: As4+Bs3+Cs2+Ds+E=0'); 

writeln(fii,’A = ',A1); 

writeln(fli/B = ',B1); 

writeln(fli,'C = ',Ci); 

writeln(fil/D = ',Di); 

writeln(fii/E = ',Ei); 

writelnC'The characteristic eqn is: As4+Bs3+Cs2+Ds+E=0'); 
writeln('A = ',Ai); 
writeln('B = 
writelnCC = ',Ci); 
writeln('D = ',Di); 
writeln('E = ',E1); 
end; C New_ChEqn } 

procedure Roots; 
var 

Ri,R2,Ii,I2,R3,I3,R4,I4 :real; 

Al,Bi,Cl,Di,Ei:real; 
procedure Real_roots(Rii:real); 
var 

t_half,t_double :real; 

begin C Real_Roots } 
if (Rii < 0.0) then 
begin 

t_half :=(0 .69*t_star)/abs(R i 1); 
writeln(fll/t_half = ',t_half :7:3/ sec'); 
writeln('t_half . = ',t_half:7:3,' sec'); 
end 
else 
begin 

t_double ;=(0 .BSjf t_star)/ abs(Rii); 
writeln(fii,'t_double = ',t_double:7:3/ sec'); 
writeln('t_double = ',t_double:7:3,' sec'); 
end; 

end; £ Real_Roots } 

procedure Find_F requency(R 1 l,R22,li 1,122 :real); 



var 


Period_T,t_half,N_half :real; 

t_double;N_clouble real; 

begin C Find_Frequency } 
if (Hi 0 0.0) then 
begin 

Feriod_T;=(2*pi*t_star)/abs(Iii); 
writeln(fii/Period_T = ',Period_T:7;3/ sec'); 
writeln('Period_T = ',Period_T :7 :3/ sec'); 
if (Rii < 0.0) then 
begin 

t_half:=(0.69*t_star)/abs(Rli); 
N_half:=t_half/Period_T; 
writeln(fii/t_half = ',t_half :7:3,' sec'); 
writeln(fii/N_half = ',N_half ;7:3,’ cycles'); 
writeln('t_half = ',t_half ;7:3/ sec'); 
writeln('N_half = ',N_half :7:3,' cycles'); 
end 
else 
begin 

t_double ;=(0 .69*t_star)/abs(Rii); 
N_double:=t_double/Period_T; 
writeln(fil,'t_double = ',t_double:7:3,' sec'); 
writeln(fil/N_double = ',N_double:7;3/ sec'); 
writeln('t_double = ',t_double:7:3,' sec'); 
writeln('N_double = ',N_double:7:3/ sec'); 
end; 
end 
else 
begin 

Real_Roots(Rii>; 

Real_Roots(R22); 

end; 

end; C Find_Frequency } 
begin C Roots ) 

New_ChEan(Ai.<Bi»Cl,DiiEi); 



Ferrari(Al,Bi,Ci,Di,El,Ri,R2,R3,R4,Il,I2,I3,I4); 

writeln(f 11/ 

writeln(fll); 

writeln(f 11/The roots of the oh eqn are 
writelnCThe roots of the oh eqn are 
writeln(fll/Lannbdal = ',R1/ + i (',11/)'); 
writeln(f ll,'Lambda2 = ',R2,' + i (',12,')'); 
writeln('Lambdal = ',R1,' + i (',11,')'); 
writeln('Lambda2 = ',R2,' + i (',12,')'); 
writeln(f 11,'For the above pair of roots :'); 
writeln('For the above pair of roots :'); 

Find_Fr equency (R i,R2,l 1,12); 
writeln(fll,'Lambda3 = ',R3,' + i (',13,')'); 
writeln(fll,'Lambda4 = ',R4,' + i (',14,')'); 
writeln('Lambda3 = ',R3,' + i (',13,')'); 
writeln('Lambda4 = ',R4,' + i (',14,')'); 
writeln(f 11,'For the above pair of roots :'); 
writeln('For the above pair of roots :'); 

F ind_Frequency(R3,R4,13,14); 
writeln(f 11); 
end; C Roots ) 

procedure Read_Gust_Data; 
begin 

writeln(output,'Input data for step gust'); 
writeln(output,'Max Gust Vel Umax = '); 
readln(input,Wmax); 

Wg :=Wmax/PlaneProg[ Ve]; 

FWg:=0; (Derivative of Gust Velocity ) 
end; 

procedure Sin_gu5t; 
begin 

writeln(output, 'Input data for sinusoidal gust'); 
writeln(output,'Max Amplitude (Gust Vel) Wmax = '); 


writeln(output, 'Input freauency(rad/s) of sinusoidal gust')- 
readln(input,F requency); 

Wg :=(Wmax/PlaneProg[ VeDf sin(F requencyf time); 
end; 

procedure RK( Y matrix; 

var F. matrix); 
begin C RK } 

FCw]:=(-(AgfYCu3+AifYCw3+AjIY[q3+Ak*Y[theta])+GvfFWg+Gw*Wg+ 

E2)lcelevator_deflection)/Ah; 

F[u3 ~(“(Ab3KYCu3+Ac^FCw3+Ad7lcYCw3+Ae)lcYCq3+Af7lfYttheta3)+GtI^FWg+ 

Gu^Wg+Ei7(celevator_def lection)/ Aa; 

FCq3:=(-(AlItY[u3+Am5KFCw3+Ao)KYCw3+Aq)lfYCq3+Ar)KY[theta3)+GxIKFWg+ 

Gy7KWg+E3IKelevator_deflection)/Ap; 

FCtheta3.-=YCq3; 
end; { RK 3 

procedure RKM(Delta_ t,elevator_def lection real; var Y .matrix); 
var 

I :perturbation_quantity; 

Y1,K1X2,K3,K4 matrix; 

T emporary .matrix; 

Dw real; 

begin C of RKM > 
for I:=u to theta do YiCI3.-=Y[I3; 

RK(Yi,Ki); 

for hssu to theta do 
begin 

Ki[I3:=delta_tIKi[13; 

Yi[I3:=YiCI3+(Kin3/2); 

end; 

RK(Y1X2); 

for I;=u to theta do 
begin 

K2[I3:=delta_tIK203; 

Yltl3:»YCI3; 



YlCI];=Yi[I]+(K2m/2); 

end; 

RK(Yi,K3); 

for I:=u to theta do 
begin 

K3[I3:=delta_tfK3[I]; 

YiCD:=Y[l]; 

YiCI3;=Yi[I]+K3CI]; 

end; 

RK(Yi,K4); 

for I:=u to theta do 
begin 

K4CI3:=delta_t*K4CD; 

YCI3:=Y[I3+(l/6)*(KiCI3+2*K2[13+2*K3[13+K4[I]); 

end; 

RKCY,T emporary); 

Dw ;=T emporaryCw3; 

Load_f actor :=-(sar(PlaneProgCVe3)*(Dw-YCq3))/(Ch_length*9r); 
writeln(fii,time,' ',YCu3,' ',Y[w3, ’ ',YCq3,' ’,Y[theta],' ',Load_Factor); 
writelnCtime/ ',YCu3/ ',Y[w3, ' ',YCq3,' ',Yttheta3/ ',Load_Factor); 
end; C of RKM } 

procedure Print_start(Y:matrix); 
begin 

writeln(fii,time,' tY[u3/ ',Y[w3, ' ',YCq3/ ’,Y[theta3,' ',Load_Factor); 
writeln(time,’ ',Y[u3,' ',YCw3, ' ',Y[q3/ '.Y[theta3/ Moad_Factor); 
end; 

procedure Print_Reqd_Data; 
begin C Print_Reqd_Data 3 
writeln(fii); 

writeln(fii,' ':iO/RESPONSE TO STEP GUSTS'); 

writeln(fii/ 'dO,' 

writeln(f li); 
writeln(fil); 

writeln(fll/An aircraft flying at an altitude of ',Initial_altitude:6:i,' ft 



writeln(fii,'at a speed of ',PlaneProgCVe3:4:i,' ft/sec encounters 
writeln(fii); 

writeln(fii,' Step Gust (ft/sec): W = ',Winax:4:2,' = constant'); 

writeln(fli); 

writeln(fll); 

writeln(fli/ ';2,'Time(sec)'/ ':6,'u_hat'/ ':8/w_hat'/ ';8, 

'a_hat'/ ':8/theta'/ ':8,'Ld_fctr'); 
writeln(' ' :2/Time(sec)',' ':6/u_hat',' ':8,'w_hatV '••8, 

'q_hat'/ ';8, 'theta',' ':8,'Ld_fctr'); 
end; I Print_Reqd_Data } 

procedure Initial_conditions; 
var 

i :perturbation_quantity; 
begin 

for l:=u to theta do Response[I3:=0.0; 
end; 

function power( x,y:real ) real; 
begin C power } 
power := exp(y*ln(x)>; 
end; C power } 

procedure Read„the_Gust; 
begin 

writeln('lnput the value of L'); 
readln(input,Scale_turbulence); 
writeln('Input the value of Sigma_w'); 
readln(input,Sigma_w); 
end; 

procedure Von_Karman; 
var 

A,B,C,D real; 
begin C Von_Karman } 

A i.339fScale_turbulence*big_omega; 


gust as be 


B ;= i+((8^fcSqr(A))/3); 

C := Power((i+Sqr(A))/il/6)); 

D ;= (Sqr(Sigma_w)*Scale_turbulence)/pi; 
psd_w ;= (D7l(B)/(C); 
end; C Von_Karman } 

procedure Et_Gust_Model; 
var 

A,B,C ;real; 

begin £ Et_Gust_Model } 
A;=:i+3fsqr(big_omegaIScale_turbulence); 

B :=sqr( i+sqr(big_onnega*Scale_turbulence)); 
C := (Sqr(Sigma_w)^(tScale_turbulence)/(2*pi); 
psd_w;= (C*A)/B; 
end; C El_Gust_Model ) 

procedure New_TFs; 
var 

Den_real,Den_imag .real; 

Num_l_real,Num_i_imag real; 
Num_2_real,Num_2_imag real; 
Num_3_real,Num_3_imag real; 
Num_LF„real,Num_LF_imag real; 

LF_coeff real; 

Num_4_real,Num_4_imag real; 
Num„5_raal,Num_5_ifnag real; 
Num_6_real,Num_6„imag real; 
Num_7_realjNum_7_imag real; 
Num_8_real,Num_8_imag real; 
Num_LFe_.real real; 

Num_LFe_.imag real; 

Common_Den ' real; 

Num_u_gust,Num_w_gust real; 
Num_a_gust,Num_theta_gust real; 
Num_LF_gust real; 

Num_u_.elev,Num_w_elev real; 



Num_q_elev,Num_theta_elev real; 

Num_LF_elev real; 

Factor real; 

begin £ New_TFs } 

£ Common denominator for the transfer functions } 

Den_real ••=(Power (F requency,4)*(Aa* Ap* Ah))+ 
(-Power(Frequency;2)*(Ab*(Ah*Aq+Ap*Ai-Am*Aj)+Aa*(Ah*Ar+Ai* 
Aq-Am*Ak-Ao«Ai)-Ag*(Ac*Aq+Ad*Ap-Am*Ae)+AllC(Ac*Aj-Ah*AB))) 
+(Ab*(AifAr-AofAk)-Ag*(Ad*Ar-AoIAf)+Alf(Ad*Ak-Ai*Af)); 
Den_imag:=(-Power(FreQuency3*(Ab*ApIAh+Aa*(Ah*Aq+ApIAi-Am*Aj)- 
Ag^tAc*Ap))+ 

(Frequency*(Ab*(AhIAr+Ai*Aq-Am*Ak-Ao*Aj)+Aa*(Ai*Ar-Ao*Ak)- 
Agf(Ac^Ar+AdilfAq-Am!fAf-Ao*Ae)+AlI(Ac*Ak+AdIAj-Ahf Af- 
Ae*Ai))); 

£ The Gust transfer functions } 

Num_ i_r eal := (power{F requency,4)I(-Gvf Ap* Ac+Gtl Ap* Ah))+ 

(-power(Frequency,2)*(Gu*(ApIAi+Ah*Aq-Am*Aj)-Gw*(AcfAq+ 

Ad*Ap-Am*Ae)-Gv*<Ac*Ar+Ad*Aq-Am*Af-Ao«Ae)+Gy*(Ac*A 3 - 

Ah«Ae)+Gx5K(Ac*Ak+AdfAj-Ah*Af-AefAi)+Gt*(Ar*Ah+AiIAq- 

Am*Ak-Ao^tAj)))+ 

(Gu*(AiIAr-Ao*Ak)-Gw^ie(Ad*Ar-Ao*Af)+Gy*(Ad*Ak-Ai*Af)); 

Num_ i_imag ;=(-power<F requency,3)*<Guf Ap^lAh-Gw*Ap* Ac-Gv^K(Acl Aq+Ad*Ap- 
Am*Ae)+Gx*(Ac*Aj-Ah*Ae)+Gt*(Ap*Ai+Ah*Aq-AmIAj)))+ 
(Frequency*(Gu*(Ar*Ah+AiIAq-Am*Ak-Ao*Aj)-Gw*(Ac*Ar+Ad*Aq 
-Am*Af-Ao*Ae)-Gv*<Ad*Ar-Ao*Af)+Gy*(Ac*Ak+Ad*A3-Ah*Af- 
Ae*Ai)+Gx*(AdIAk-Ai*Af)+Gt*(Ai*Ar-Ao*Ak))); 

Num_2„real :=(power(F requency,4)*(Aa*Gv*Ap))+ 

(-power(Frequency,2)^(c(Aa5K(Gv*Ar+Gwf Aq-Gxf Ak-Gy*Aj)+ 
Abf(GwfAp+Gv*Aq-GxfAj)-Ag*(GufAp-Gx^eAe+Gt*Aq)-(Gv*Ae- 

Gt«A3)*Al))+ 

(Ab5l<Gw*Ar-GyfAk)-Ag*(Ar*Gu-Gy*Af)+Al*(Gu*Ak-GwfAf)); 

Num_2_imag;=(-POwer(Frequency,3)*<Ab*Gv*Ap+Aa*(Gw*Ap+Gv*Aq-Gx*Aj>+ 

Gl*Ap))+ 

(Freauency*(Aa*(GwfAr-Gy*Ak)+Ab*(Gv*Ar+Gw*Aq-Gx*Ak-Gy^ 

A3)-Ag«(Gu*Aq-Gx*Af-Gy*Ae+Gt*Ar)+AlI(Gu*Aj-Gv*Af-Gw*Ae 


+Gt*Ak))); 



(Numerator for Theta} 

Num_3_real;=(-power(Frequency,2)^fc(Aaf(Ah*Gy+Ai*Gx-Am*Gw-Ao*Gv)+ 
Ab*(Ah*Gx-Am*Gv)-Ag*(AcfGx-Am*Gt)+Alf(Ac*Gv-Ah*Gt)))+ 
(Ab*(Ai*Gy-Ao*Gw)-Agf(Ad*Gy-Ao*Gu)+AH(AclfGw-Ai*Gu))j 
Num_3_imag :=:(-pawer(F requency,3)^«(Aa*(Ah*Gx-Am*Gv)))+ 
(FreQuency*(Aa*(Ai*Gy-Ao*Gw)+Ab*(Ah*Gy+Ai*Gx-Am*Gw- 
Ao^Gv)-Ag*(AcfGy+Acl*Gx-Am*Gu-Ao*Gt)+Alf(Ac*Gw+Ad*Gv- 

Ah*Gu-Ai^KGt))); 

£ The elevator transfer functions) 

Num_4_real;= (-power(Frequency,2)*(Ei*(AifAp+Ah*Aq-Am*Aj)-E2I(AdiKAp+ 
AclAq-Am^lcAe)+E3*(Ac*Aj-Ah*Ae)))+ 
(Ei*(Ai*Ar-Ao*Ak)-E2*(Ad^(tAr-Ao*Af)+E3*(Ad*Ak-Ai*Af)); 

Num_4_imag := <-power(Frequency,3)*<Ap3K(Ei*Ah-E2*Ac)))+ 

(Frequency*(Ei*(Ar*Ah+Ai*Aq-Am*Ak-Ao*Ao)-E2f(AdfAq+ 

Ac*Ar-Am5tAf-Ao*Ae)+E3I(Ac*Ak+AdfAj-AhfAf-Ai*Ae))); 

Num_5_real -■= (-power(F requency,2)*(-Ei^Apf Ag+E2*<Abf Ap+Aaf Aq)- 
E3*Aa*Aj))+ 

(-Ei*(Ag*Ar-AlfAk)+E2*(Ab*Ar-Al*Af)-E3*(Ab*Ak-Ag*Af)); 

Num_5_imag •'= (-power (F requency>3)*(E2f Aaf Ap))+ 

(F^equency^K(-Eif(Agi^Aq-Alf Aj)+E2*(Aa*Ar+Ab*Aq-Al*Ae)- 
E3*(AaacAk+Ab*Ak-Ag*Ae))); 

C Numerator for ElevatorTheta ) 

Num_6_real := (-power(Frequency,2)^K(E3*Aa5lAh-E2TlcAa^((Am))+ 
(Ei5l(Ao*Ag-Al*Ai)-E2*(Ab*Ao-Al*Ad)+E3*(Ab*Ai-Ag*Ad)); 

Num_6_imag ;= (F requencyit(Eii(((Am^K Ag-AlilcAh)-E2^lf (Aa^ltAo+Ab^f Am-Ac5lcAl)+ 

E3*(Ab*Ah+Aa«Ai-Ac*Ag))); 

( Numerator of TF q (gust) ) 

Num_7_real :» -F requency f Num_3_imag; 

Num_7_imag ;= F requencyf Num_3_real; 

£ Numerator of TF q (elev) ) 

Num_8_real -F requency^cNum_6_imag; 

Num_8_imag := F requency^KNum_6_real; 

LF_coeff := -(2*sqr(PlaneProgCVe3))/(gr*PlaneProg[c3); 

Num_LF_real:= LF_cosff«((-Frequency*Num_2_imag)-Num_7_real); 
Num_LF_imag := LF_coef f *((Frequency*Num_2_real)-Num_7_imag), 

Num_LFe_real.« (.gr/PlaneProg[CLe3)*((-Frequency)ieNum_5_imag)-Num_8_real), 


Num_LFe_imag:= <-9r/PlanePrD9[CLsB*((Freciuencg*Num.5_reaU-Num_8_imag), 

C Modify Common denominator for Transfer Function theta 3 
Comfnon_Den := sqr(Den_real)+sqr(Den_imag); 

Num_u_gust := sqr(Nurn_ i_real)+sqr(Num_ i_imag); 

Num_w_gust := sqr(Num_2_real)+sqr(Num_2_imag); 

Num_q_gust := sqr(Num_7_real)+sqr(Num_7_imag); 

Num_theta_gust:= sqr<Num_3_real)+sqr(Num_3_imag); 

Num_LF_gust := sqr(Nurn_LF_real)+sqr(Num_LF_imag); 

Num_u_elev ;= sqr(Num_4_real)+sqr(Hum_4_imag); 

Num_w_elev ;= sqr(Num_5_real)+sqr(Num_5_imag); 

Num_q_elev ;= sqr(Num„8_real)+sqr(Num_8_imag); 

Num_theta_elev := sqr(Num_6_real)+sQr(Num_6_imag); 

Num_LF_elev ;= 5qr(Num_LFe_real)+sqr(Num_LFe_imag); 

Factor := i/PlaneProg[Ve3; 

TF_u_gusi ;= F actor^csqrt(Num_u_gust/Common_Den); 

TF_w_gust ;= F actorjif sqrt(Num_w_gust/Common_Den); 

TF_q_gust := Factor* sqrl(Num_q_gust/Common_Den); 

TF_theta_gust := F actor*sqrt(Num_theta_gust/Common_Den)j 
TF_LF_gust ;= Factor*sqrt(Nunn_LF_gust/Common_Den); 

Ph„u„gust:= (i80/pi)*<arctan(Num_i_imag/Num_l_real)- 
arctan<Den_imag/Den_real)); 

Ph_w_gust :o (180/pi)*(arotan(Num_2_imag/Num_2_real)- 
arctan(Den_iinag/Den_real)); 

Ph_q.gust ••=« (i80/pi)f (arctan(Num_7_imag/Nuni_7_real)- 
arctan<D«n_imag/Den_real)); 

Ph_thela_gust ;= (180/pi)*(arctan(Num_3_imag/Num_3_real)- 
arctanCDen>.imag/Den_real)); 

Ph_LF_gust:= <i80/pi)*<arctan(Num_LF_imag/Num_LF_real)- 
arctan(Den_imag/Den_real)); 

TF_u_elav sqrt(Nuni_u_elev/Common_Den); 

TF_ w_elev := sqr t(Num_w_elev/Common_Den); 

TF„q_elev ;=: sqrt(Num_q_ele v/Common_Den); 

TF_theta_elev:= sqrt<Num_theta_elev/Common_Den); 

TF_LF_elev := sqr t(Num_LF_elev/ Common_Den); 

Ph_u_elc V := (l80/pi)*<arctan(Num_4_imag/Nuin_4_real)- 
arotan<Den_imag/Den_real)); 


Ph_w_elev:= (i80/pi)*(arctan(Num_5_ima9/Num_5_real)- 
arctan(Den_imag/Den_real)); 

Ph_q_elev ;= (180/pi)f (arctan(Num_8_imag/Num_8_real)- 
arctan(Den_imag/Den_real)); 

Ph_theta_elev := (i80/pi)5t(arctan(Num_6_imag/Num_6_real)- 
arctan(Den_imag/Den_real)); 

Ph_LF_elev := (i80/pi)*(arctan(Num_LF e_imag/Num_LFe_real)- 
arctan(Den_imag/Den_real)); 

TF_u_gust:= sqr(TF_u_gust); 

TF_w_gust:= sqr(TF_w_gust); 

TF_q_gust:= sqr(TF_q_gust); 

TF_theta_gust ;= sqr(TF_theta_gust); 

TF_LF_gusl;= sqr(TF_LF_gust); 
end; £ New_TFs 3 

procedure Get_Gain_Factors (var pr:glnp); 
begin C Get_Gain_Factors 3 
read(kil,prEi3); 
read(i<ii,pr[23); 
read(kii,pr[33); 
read(kiijpr[43); 
read<kli,prC53); 
read(kii,prC63); 
read<kli,prC73); 
readln<kil); 

writeln(fii,'Kalpha_gust = ',prCi3); 
wrileln(flVT_Filter = ',pr[2]); 
writeln(fll,'Ku = ',prC33); 
writeln(f ii,'Kalpha = ',prC43); 
wrileln(f ll/Kq = ',prC53); 
writeln(fll,'K theta = ',pr[63); 
writelnCfli/KLdFctr = ',prC73); 
end; £ Get_Gain_Factors 3 

procedure Gain_F actors(pr -’glnp); 
begin 


Kalpha_gust :=pr[i]; 

T_Filter;=pr[2]; 

K;u;=pr[33; 

Kalpha:=pr[4]; 

Kq;=prC5]; 

Ktheta ;=prC6]; 

KLdFctr:=pr[7]; 

end; 

procedure Out_spectra; 
begin 

Out_u_gust ~TF_u_gustf psd_w; 

Out_w_gust :=TF_ w_gust*osd_w; 

Oul_q_gust :=TF_q_gust5lcp5d_w; 
Out_theta_gust:=TF_theta_gustIpsd_w; 

Out_LF_gust ■•=TF_LF_gustiKpsd_w; 
end; 

procedure Print_ing; 
begin 

writeln(fll,' ','L = ',Scale_turbulence:3:2/ ':5, 'Sigma = 
Sigma„w:3:3); 

writelnCfii/ ':2,'big_otinega',' '.•6/Frequency'/ ':3/PSD_W', 

' ':9/TFu'/ ';iO/TFw'/ ';iO/TFq'/ '.-V/TF.LdFctr', 

' ':6/0ut^u'/ ':8/0ut.w'/ ':8/0ut_q'/ ':6, 

'Out„LdFclr'); 

writeln(fll/ ':3/(rad/ft)'/ '.'S/Wz)'); 

writeln(pil/THE TRANSFER FUNCTIONS FDR THE ELEVATOR'); 
writeln<pii/ '/L = ', Scale. turbulence :3:2/ 'S/Sigma = 
Sigma_.w:3:3); 

writelnCpll/ ':2/big_omega'/ ':6/Frequency'/ ';3/PSD_W', 

' ':9/TFuV ':iO/TFw'/ ':10/TFtheta'/ '=3); 
writelnCpii/ ':3/(rad/ft)'/ ';6/(Hz)'>; 
writelnC '/L = ',Scale_.turbulence:3;2/ ':5/Sigma = 
Sigma.w:3:3); 

writelnC' ';2/big_omega'/ ':6/Frequency'/ ':3, PSD_W, 



' ':9/TFuV ':iO;TFwV MO/TFq',' '--T/TFadFctr', 

' ';6/0ut_uV ':8/Out_w'/ ':8,'0ut_qV ’:6, 

'Out_LdFctr'); 

writelnC' ':3,'(rad/fl)',' '.-6, '(Hz)'); 
end; 

procedure Print_out; 
begin 

writeln(fil,big_omega:iO,' ',dim_frequency:10/ ',psd_w;iO/ 
TF_u_gust:10,' ',TF_w_gust;iO,' 'JF_q_gust:iO,' 
TF_LF_gusl:iO/ ',0ut_u_gust;10/ ',Out_w_gust:iO/ ', 
Out_q_gust ; 10/ ',Out_LF_gust -iO); 

writeln(pll,big_onnega:iO/ ',dim_frequency;iO/ ',psd_w:iO/ ', 
TF_u_elev:iO/ ',TF_w_elev:10/ ',TF_theta_elev:iO); 
writeln(big_omega:10/ ',dim_frequency;10/ ',psd_w:10/ 
TF_u_gust:iO/ ';TF_w_gust:iO/ 'JF_q_gust:iO/ 
TF_LF_gust:iO/ ',Out_u_gust:iO/ ',Out_w_gust:iO/ ', 
Out_q_gust : 10/ ',Out_LF_gust :10); 
end; 

procedure Graf_Print_out; 
begin 

writeln(kkl,big_ome 9 a/ ',psd_w); 

writeln(l<k2,bi9_omega/ ',TF_u_gust/ ',TF_w_gust/ ',TF_q_ 9 ust, 
' MF_LF_gust); 

writeln(kk3,bi9_ome9a/ ',0ut_u_9ust/ ’,Out_w_gust/ ', 
Oub_q_gust/ ',0ut_LF_9ust); 

wrileln(kk4,bi9_omega/ ',Frequency/ ',TF_u_elev/ ',TF_w_elev, 

' ',TF_theta_elev); 
end; 

procedure IntegrateCik ■•integer;Out_quantity ;real; 

var sum:real); 
begin C Integrate } 
if ((ik = 1) or (ik = 899)) then 
sura ;=sum+Out_ quantity 


else 

if ((ik mod 2) = 0) then 
sum :=sum+4^IOut_quantity 
else 

sum;=surn+2f0ut_quantity; 
end; C Integrate 3 

procedure Relax_St_Stblty; 
begin { Relax_Static_Stability 3 
CwritelnCGive the value of KFactor_CLaw'); 
readln(KFactor_CLaw);3 
KF actor_CLaw :=□; 

ACMtheta :=KF actor_CLaw*PlaneProg[CMalpha]; 

ACMalpha :=(i-KF actor_CLaw)IPlaneProg[CMalpha]; 
writeln(fii); 

writeln(fii,' ') 

writeln(fii/KFactor_CLaw ='XFactor_CLaw); 
end; (Relax_Static_Stabilty 3 

procedure Step_gust_calc; 
begin C Step_gust_calculations 3 
Read_Gust_Data; 

Print_Reqd_Data; 

Initial^conditions; 

time:=0.0; 

Load„F actor ■•sO .0; 

Print_start(Response); 
time :=timB+delta_ time; 
repeat 

RKM«delta_time/t_star),elevator_deflection, response); 

time :=time+delta_tiroe; 
until (time > Final_time); 
end; C Step_gust_calculations 3 


procedure Check„Print(ji integer); 
begin C Check_pr'irtt 3 



if ((oi=i) or (ji=10) or (ji=100)) then 
count :=i00000/ji; 
check_coef f ~big_omegaicount; 

if (((check_coeff-trunc(check_coeff)) = 0.0) or (big_omega=3e-4) 
or (big_ornega=6e-4) or (big_omega=7e-4)) then 
begin 
Print_out; 

Graf_Print_out; 

end; 

end; C Check_Print ] 

function func(pr;glnp);real; 
var 

jj integer; 
begin £ func } 

Gain_F actors(pr ); 

Dummy_Coeffs; 

sum2:=0; 

delta_omega:=PlaneProgCVe]/i00000; 
for jj := 1 to 999 do 
begin 

big_omega := jj/iOOOOO; 

dini_freQuency;=big„omega*PlaneProg[Ve]; 

Von_Karman; 

£Et_Gust_Model;) 

Frequency :=big_omegaf Ch_length; 

New„TFs; 

□ut_SDectra; 

lntegrate(jij0ut_LF„gust,sum2>; 

end; 

cornfort_criterion :=:sqrt((delta_omegaitsu[ii2)/(3)Kpi))j 
f unc :=comf ort_criterion; 
end; £ func } 

PROCEDURE amoeba(VAR p: glmpnp; VAR y: gimp; ndim: integer; 
ftol: real; VAR iter; integer); 


(f Programs using routine AMOEBA must supply an external function 
func(pr;glnp):real whose minimum is to be found. They must 
also define types 
TYPE 

glmpnp = ARRAY [i..mp,i..np] OF real; 
gimp = ARRAY [i..mp] OF real; 
glnp = ARRAY [i..np3 OF real; 
where mp and np are physical dimensions *) 
label 990; 

CONST 

alpha=1.0; 

beta=0.5; 

gamma=2.0; 

itmax=500; 

VAR 

mpts,j,inhi,ilo,ihi,i: integer; 
yprr,ypr,rtol: real; 
pr,prr,pbar: glnp; 

BEGIN { Amoeba ) 
writelnCEntering AMOEBA'); 
mpts ■■= ndim+i; 
iter ;= 0; 

WHILE true DO 
BEGIN 
ilo := i; 

IF (yCil > yC23) THEN 
BEGIN 
ihi := i; 
inhi ••= 2 
END 
ELSE 
BEGIN 
ihi ;= 2; 
inhi := 1 
END; 

FOR i := i to mpts DO 


BEGIN 

F (yCil < yCiloD) THEN ilo := i; 

IF (y[i3 > yCihil) THEN 
BEGIN 
inhi ;= ihi; 
ihi ;= i 
END 
ELSE 

IF (y[i3 > y[inhi3) THEN 
IF (i 0 ihi) THEN inhi := i 
END; 

rtol := 2.0fabs(y[ihi3-y[ilo3)/(ab5(y[ihi3)+abs(yCilo3)); 
IF (rtol < ftol) THEN GOTO 990; 

IF (iter = itmax) THEN 
BEGIN 

writeln('pause in AMOEBA - too many iterations'); 
goto 990; 

END; 

iter ;= iter+1; 

FOR j ;= i to ndim DO pbarCjl ;= 0.0; 

FOR i :=: i to mpts DO 
IF (i 0 ihi) THEN 
FOR j := 1 to ndim DO 
pbar[j3 ;= pbarCj3+p[i,j3; 

FOR j :s i to ndim DO 
BEGIN 

pbarCjl := pbarCj3/ndim; 
pr[j3 := (i.O+alpha)*pbar[j3-alphaIKp[ihi,33 
END; 

ypr := func(pr); 

IF (ypr <= yCilol) THEN 
BEGIN 

FOR j := i to ndim DO 

prr[j3 ~ gammaHEprCjl+d .0-gamma)iKpbar[j]j 

yprr ■•= func(prr); 

IF (yprr < yCilol) THEN 


BEGIh4 

FOR j := i to ndim DO pEihi.j] := prKj]; 
yCihi] := yprr 
END 
ELSE 
BEGIN 

FOR j := i to ndim DO p[ihi,j] := prCj]; 
y[ihi3 := ypr 
END 
END 
ELSE 

IF (ypr >= yCinhil) THEN 
BEGIN 

IF (ypr < y[ihi]) THEN 
BEGIN 

FOR 0 ;= 1 to ndim DO p[ihi,j3 := prCA 
yCihi3 := ypr 
END; 

FOR j •■= I to ndim DO prrCj3 ;= beta*pCihi,33+(i.O-beta)*pbar[33; 
yprr := func(prr); 

IF (yprr < yCihi3) THEN 
BEGIN 

FOR j := i to ndim DO p[ihi,j3 ~ prrCj3; 
y[ihi3 := yprr 
END 
ELSE 
BEGIN 

FOR i := i to mpts DO 
IF (i 0 ilo) THEN 

BEGIN 

FOR j := 1 to ndim DO 

BEGIN 

prCj3 ~ O.53t(p[i)33+p[ilo,j3); 
p[i43 := prtjl 
END; 

y[i3 := func(pr) 



END 


END 

END 

ELSE 

BEGIN 

FOR i := i to ndim DO plihi.j] := pKj]; 
yCihi] := ypr 
END 
END; 

390;END; C Amoeba 3 

FUNCTION ranKVAR idum,glixi,glix2,glix3 integer; 
var glr;galr): real; 

(tK Programs using RANi must declare the -following variables 
VAR 

glixi,glix2,glix3: integer; 
glr: ARRAY [i..973 OF real; 
in the main program. *) 

CONST 

m 1=259200; 

ial=7141; 

ici=54773; 

rmi=3.8580247e-6; (* i.O/mi *) 

m2=134456; 

ia2=812i; 

ic2=284ii; 

rm2=7.4373773e-6; i.0/m2 *) 

m3=*243000; 

ia3=456i; 

ic3=5i349; 

VAR 

jjm: integer; 

BEGIN 

IF (idum < 0) THEN 
BEGIN 

glixi •■= (ici-idum) MOD ml; 


glixi ■•= (iai*glixi+ici) HOD mi; 
glix2 := glixi MOD m2; 
glixi ;= (iaifglixi+ici) MOD mi; 
glixS := glixi MOD m3; 

FOR jjm := i to 97 DO 
BEGIN 

glixi ;= (iaiitglixi+ici) MOD mi; 
glix2 ;= (ia2*glix2+ic2) MOD m2; 
glrCjjm] := (glixi+glix2*rm2)*rmi 
END; 

idum := i 
END; 
repeat 

glixi := (ialfglixi+ici) MOD mi; 

glix2 ;= (ia2*glix2+ic2) MOD m2; 

glix3 := (ia3*glix3+ic3) MOD m3; 

jjm ■•= i + (97f glix3) DIV m3; 

until ((jjm < 97) AND (jjm > i)); 
rani ;= glrCjjml; 

glrCjjm] (glixi+glix2*rm2)Irmi; 

END;£ of RANi 3 

PROCEDURE memcof(data: glnarray; n,m: integer; VAR pm: real; 

VAR cof : glmarray; wki,wk2: glnarray; wkm; glmarray); 

(Kt Programs using routine MEMCOF must define the data types 
TYPE 

glnarray « ARRAY [i..n] OF real; 
glmarray = ARRAY Ci..m3 OF real; 

and must provide workspace arrays wki,wk2,wkm with the dimensions 
shown in the arguments above, f) 

LABEL 991; 

VAR 

k,j4: integer; 
pneum,p,denom: real; 

BEGIN 


for 0 := i to n do 

BEGIN 

p p+SQr'Icis.t/St 

END; 

pm := p/n; 
wkilil := dataCi]; 
wk2[n-i] := dataCn]; 

for j := 2 to n-i DO 
BEGIN 

wkiLi] := datalj]; 
wk2Cj-i3 := dataCj] 

END; 

FOR k ;= i to m DO 
BEGIN 

pneum — 0.0; 
denom := 0.0; 

FOR 0 := i to n-k DO 
BEGIN 

pneum ~ pnsum+wkiCo]Iltwk2tj3; 
denom := denom+sar(wkiCj])+5dr'(wk2[3]) 

END; 

cofCk3 := Z.Ofpneum/denom; 
pm := pm2(c(i.O-sqr(cof[k3)); 

IF (k 0 i) THEN 
BEGIN 

FOR i := i to k-i DO 
BEGIN 

cofCi3 := wkm[i3'00ftk3*wkmCk-i3 

END 

END; 

IF (k = m) THEN GOTO 99i; 

FOR i := i to k DO 
BEGIN 

wkmCi3 := cofCi3 
END; 

FOR j := i to n-k-i DO 



BEGIN 

wkiCj] := wkiCj3-wkm[kHwk2[o]; 
wkZCj] := wk2[j+i]-wkmCk]fwkiCj+i3 
END 
ENDj 

99i:END; C of MEMCOF > 

FUNCTION evlmem(fdt: real; cof ; glmarray; m: integer; pm: real); real 

GK Programs using routine EVLMEM must define the types 

TYPE 

glmarray = ARRAY Ci..m] OF real; 
where m is the dimension of the array of coefficients. *) 

VAR 

wr,wi,wpr,wpi,wtemp, theta : double; 
sumi,sumr: real; 
i: integer; 

BEGIN 

theta := 6.283i85307i7959*fdt; 

wpr ;= cosCtheta); 

wpi := sin(theta); 

wr := i.O; 

wi ;= 0.0; 

sumr := i.O; 

sumi ;= 0.0; 

FOR i ;= i to m DO BEGIN 
wtemp ■•= wr; 
wr ;= wrlfwpr-wiHfwpi; 
wi := wi*wpr+wtemp*wpi; 
sumr := sumr-cof[i3*sngl(wr); 
sumi := sumi-cof[i3*sngl(wi) 

END; 

evlmem := pm/(sqr(sumr)+sqr(sumi)) 

END; £ of EVLMEM 3 

procedure Integrate_psd_w; 


var 



kkj .'integer; 

begin C of Integrate_psd_w } 

5umi:=0; 

delta_omega:=PlaneProgCVe]/iOOOOO; 
for kkj:=i to 999 do 
begin 

big_omega :=kki/100000; 

Von_Karman; 

{Et_Gust_Model;) 

Integrate(kkj,psd_w,sumi); 

end; 

sum! :=suml*delta_omega; 
end; C of Integrate_psd_w } 

procedure Rand_Phi; 
begin C of Random_Phi } 

random_Phi;=2Ipi*RANi(lphi,gplixi,gplix2,gplix3,gphi); 
end; C of Random_Phi } 

procedure Generate_random(var random_omega real; var accept ;boolean); 
var 

Compare_f unCjT otal_area_A --real; 

random_Aprime,random_y real; 

begin C Generate_random 3 
Compare_f unc -'=3 .54336E4/ sumi; 

Total_area_A :=Compare_f unc*(F inal_Qmega-Initial_omega); 

random_Aprime:=Total_area_A^SRANi(Iomega,golixi,golix2,golix3,gomega); 

randorn_ornega ;=(random_Aprime/ Co(npare_f unc)+Initial_ornega, 
random_y:=RANlCIrandom_y,gylixl,gylix2,gylixagrandom_y); 

big_omega ;=random_omega; 

Von_Karman; 
psd_w :=psd_iM/sufni; 

if ((psd_w/Compare_func) <= random_y) then 
accept :=true; 

random_ofriega:=randorn_omegafPlaneProgCVe3; 

end; C Generate_randorn 3 



var 


surnrri :realj 
kkrn unteger; 
begin { Rand_Time_Func 3 
5urrim.=0; 

for kkrn—i to Total_N do 
surnm:=summ+cos(Omeg3kCkkm3ftime+PhikCkkm3); 
randon-i_T ime_F unc ;=sqrt(sumilZ'/T otal_N)fsunrim; 
end; t Rand_Time_Func 3 

PROCEDURE fourKVAR data: gldarray; nn,isign; integer); 

Programs using routine FOUR! must define type 
TYPE 

gldarray = ARRAY [i..nn2] OF real; 
in the calling routine, where nn2=nn+nn. 

VAR 

iijj j,n,mmax,m, j,istep,i : integer; 
wtemp,wr,wpr,wpi,wi, theta : double; 
tempr, tempi: real; 

BEGIN 

n := 2f nn; 

} := i; 

FOR ii := i to nn DO BEGIN 
i := 2I^ii-l; 

IF (j > i) THEN BEGIN 
tempr := dataCjl; 
tempi := dataCj+il; 
dataEj] := dataCil, 
data[j+il := dataCi+il; 
dataCil := tempr; 
dataCi+il ■= tempi 
END; 

rn := n DIV 2; 

WHILE ((m >= 2) AND (j > m)) DO BEGIN 
j := 3-rn; 


m := m DIV 2 



END; 


j ■•= j+rn 
END, 

mmax := Z, 

WHILE (n ) rnmax) DO BEGIN 
istep := Z^rnmax; 

theta .= 6.233iS530717959/(isign*mmax); 
wpr := -2.0*sqr(5in(;0.5*theta)); 
wpi := sin(theta); 
wr .= i.O; 
wi ;= 0.0; 

FOR ii ;=: i to (mmax DIV 2) DO BEGIN 
m ■■= 2fii-i; 

FOR jj ;= 0 to ((n-m) DIV istep) DO BEGIN 
i := m + jj^istep; 
j .= i+rnmax; 

ternpr •= sngl(wr)fdataCj]-sngl(wi)*data[j+l]; 
tempi ■•= sngl(wr)*dataCj+13+sngl(wi)*dataCj]; 
dataCj] := dataCil-ternpr; 
dataCj+i] := dataCi+i]-tempi; 
dataCi] := dataCi]+tempr; 
dataCi+i] := dataCi+il+tempi 
END; 

wtemp ;= wr; 

wr ;= wr*wpr-wif wpi+wr; 
wi := wi*wpr+wtemp^fewpi+wi 
END; 

mmax ;= istep 
END 
END; 

begin { main } 

Reset(DskFl/AirplaneData'); 

{Reset(DskFl/AirplaneEtkin');} 

Re5et(inp/vertex .in'); 

ResetCkli/gain.in'); 



Rswrite(f ii/tfg.o'); 

Rewrite(pi i/tf e .o'); 

Rewrite(kkl,'gpsdw .dat'); 

Rewrite(kk2/gtf g .dat'); 

Rewrite(kk3/goutg .dat'); 

Rewrite(kk4,'gtf e .dat'); 

RewriteCkk5/'gtime.dat'); 

ReadDskForACDer; 

Am_ConstCoeff; 

Read_the_Gust; C for Von Karman model of gust ) 

driver for OPTIMISATION 

C writelnC'Reading the vertices of the initial simplex'); 
for i:= i to mp do 
begin 

for j;=l to np do 
begin 

read(inp,p[i,j])i 

end; 

readln(inp); 

end; 

ndim ;= np; 

FOR i := i to mp DO 
BEGIN 

FOR j := i to np DO xCjl -'= pCiJl; 
yCi] := func(x); 

END; 

amoeba<p,y,ndim,ftolAter); 

writeln; 

writeln('Iterations: 'jiterS); 
writeln('Vertices of final simplex and'); 
writelnt'f unction values at the vertices:'); 
writeln; 

writelni'i' :3/Kalphagust' :i3/T_F liter' :iO/Ku' :i2, 

'Kalpha':i2,'Kq':i2,'Ktheta':i2,'KLdFctr';i2, 

'function' :i4); 



wr ileln; 

FOR i := i to mp DO 
BEGIN 
i.Mrite(i G); 

FOR j := i to np DO 
BEGIN 

write(pCi,o]:i2.6); 

END; 

writeln(y[i].i2:6); 

END; 

writeln;} 

C driver tor CALCULATIONS USING CONTROL LAW 

{Relax_5t_Stblty; 

Get_Gain_F actors(pr); 

Gain_F actors(pr); 

Dummy_Coeffs; 

Roots; 

Step_gu3t_calG; 
sum2 :=0; 

delta_o(nega:=PlaneProg[Ve3/i00000; 

Print_ing; 

for j ;= i to 999 do 
begin 

big_ornega := j/iOOOOO; 

dim_frequency :=big_omegalPlaneProgCVe]; 

Von_Karrnan; 

Et_Gust_Model; 

F requency ;=big_omega*Ch_length; 

New_TFs; 

Out_Spectra; 

Check_Print(j); 

Integrate(j,Out_LF_gust,sutn2); 

end; 

comfort_criterion:=sqrt((delta_omega*sum2)/(3*pi)); 
writeln(fii/Comtort_Criterion = ',comfort_cr iter ion); 



writeln('Comfort_CriteriQn = 'jComfort_criterion);} 

<■ driver for RANDOM GUST CALCULATIONS } 

Iomega :=- 123; 

Iphi:=-il3; 

Irandom_y:=-il; 

3ccept;=false, 
lntegrate_psd_w; 
for nkki:=i to TotalJM do 
begin 
Rand_Phi; 

PhikCnkki];=random_phi; 

repeat 

Generate_random(random_Qmega, accept); 
until (accept = true), 

0rnegak[nkki3:=random_omega; 

end; 

time :=Initial_tirne; 
iicount :=0; 
repeat 

Rand_T ime_F unc; 
iicount :=iicount+i; 

T ime_dataCiicount3 ;=r andom_T ime_F unc; 
if (iicount < 200) then 
begin 

writeln(kk5,tirne,Tirne_data[iicount]); 

end; 

time :=time+delta_time; 
until ( iicount = Total_N ); 
for kkm:=i to Total_N do 
begin 

WKi[kkm3:=0; 

WK2i:kkm]:=0; 

end; 

for kkm:=i to 0rder_m do WKM[kkm3-=0; 

MEMCOF(Time_data,Total_N, Order _M,PM, COP, WKi,WK2,WKM); 



Ralax-_5t_Stblty; 

Get_Gain_Factors(pr,'; 

Gain_F actors(pr); 

Dummy_Coeffs; 
for kkm;=i to nn2 do 
begin 

Moutu[kkm]:=Q, 

MoutwCkkm] :=0, 

MoutqCkkmD .=0; 

MoutLFCkkm].=0; 

end; 

for kkm;=i to nnby2 do 
begin 

F requencys :=kkm/(nn^delta_time); 
big_omega ■=Fr equencys/Planepr ogCVe]; 
dim_f requency :=2*pi*Frequency3; 
f dt :=F requency5*delta_time; 
psd_w;=EVLMEH(fdt,COF,Order_M,PM); 
if (kkm < nnby2) then 
begin 

writelnCkki,big_omega,' ',psd_w); 
writeln(big_omega,' ',psd_w); 
end; 

F requency ;=big_omega^SCh_length; 
New_TF5; 

□ut_Spectra; 

MoutuC2Skknn+i] ;=0ut_u_gust; 
MoutwC2^l^kkm+ ii :=0ut_ w_gu3t; 

MoutqC2:f kkrin + LI ~0ut_q_gust; 
MoutLFC2*kkm+ i] :=Out_LF_gust; 
end; 

f our i (Moutu .nnisign); 
fouri(MoutWjnn,isign); 
f our iCMoutqjnn,isign>; 
fouri(MoutLF,nn,isign); 
time:=Initial_time; 



for kkrrr=i to 133 do 
bEgin 

writeln(kk3,time/ ',MouiuC2Ikkm-iV ',Moutw[2*kkm-i], 
' ',Moutq[2*kkm-i]/ ',MoutLFC2*kkm-i3); 
writelnCtime/ ',MoutuC2*kkm-i]/ ',Moutw[2Ikkm-i], 

' ',Moutq[2if:kkm-i3/ ',MoulLFC2:Skkrn-i]); 
tirns:=tirne+delta_time, 
endj 


Close(fii); 
Close(inp); 
Close(pii); 
Close(kii); 
Close(kki); 
Close(kk2); 
Close(kk3); 
Close(kk4); 
Close(kk5); 
end . { main } 
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